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Chapter  1.  Introduction 


1.1.  Motivation 

The  design  of  high  current  plasma  devices  such  as  arcjet  thrusters,  arc  heaters, 
magnetohydrodynamic  generators,  and  plasma  opening  switches  has  in  the  past  been 
based  mostly  on  empirical  methods.  In  the  case  of  arcjet  thrusters,  these  methodologies 
have  produced  a  steady,  reliable  device  when  operated  at  low  power  levels.  However,  it 
is  desirable  to  increase  the  performance  and  operating  range  of  the  arcjet  so  that  they  may 
run  more  efficiently  and  reliably  and  be  applicable  to  a  wider  set  of  missions  [1,2]. 
Because  of  the  expense  and  time  cost  of  empirical  methods,  attention  has  been  turned  to 
developing  accurate  numerical  simulations  which  aspire  to  aid  in  the  development  of 
higher  power  designs,  while  providing  a  clearer  understanding  of  internal  physical 
processes.  However,  due  to  limitations  in  computing  power  and  the  complexity  of 
plasma  behavior,  simplifications  must  be  made  to  keep  the  calculations  tractable. 

It  is  well  known  that  plasmas  are  an  inherently  unstable  medium.  Many  steady 
state  situations  are  precarious  balances  between  competing  mechanisms  which  can 
become  unstable  if  those  balances  are  slightly  upset.  The  length  and  time  scales  over 
which  these  instabilities  occur  can  be  very  small  and  difficult  to  capture  in  a  larger  scale 
simulation.  Hence,  much  of  the  important  smaller  scale  physics  may  not  be  captured  by  a 
computational  model,  resulting  in  an  inability  to  capture  real  effects.  In  addition,  the 
desire  to  calculate  steady  state  solutions  removes  any  ability  to  capture  realistic  time 
dependent  phenomena.  Existing  codes  [3-6]  for  arcjet  simulation  also  impose  a  2D, 
axisymmetric  condition.  These  codes,  while  ambitiously  trying  to  model  the  overall 
behavior  of  the  arcjet,  may  not  have  the  spatial  or  time  dependent  capabilities  to  capture 
small  scale  and  asymmetrical  effects.  As  a  result,  care  must  be  exercised  when 
interpreting  the  computed  results  for  cases  where  there  is  no  existing  experimental  data  to 
compare  with. 
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One  phenomena  that  is  sometimes  seen  during  arcjet  operation  is  the  transition 
from  a  stable,  diffuse  mode  of  arc  attachment  on  the  anode  to  a  destructive,  spotty,  and 
constricted  attachment  mode  [7,8].  This  phenomena  is  essentially  an  important  limiting 
factor  to  the  applicability  of  arcjets  to  more  demanding  space  missions.  A  schematic  of 
the  arcjet  running  in  the  diffuse  mode  is  given  in  Fig.  1.1.  A  possible  cause  of  this 
transition  is  an  instability  which  occurs  in  the  near  electrode  region,  much  like  the  way 
instabilities  in  a  laminar  flow  boundary  layer  lead  to  the  transition  to  turbulence.  These 
instabilities  would  in  theory  be  initiated  by  small  perturbations  in  plasma  properties 
which  grow  in  time.  The  purpose  of  this  study  is  therefore  to  model  the  time  dependent 
response  of  small  oscillating  disturbances  in  the  near-anode  region  and  to  identify 
mechanisms  which  affect  the  growth  or  decay  of  these  disturbances. 


Figure  1.1.  Cross  section  schematic  of  an  arcjet  running  in  the  diffuse  attachment  mode. 

The  region  near  an  anode  in  contact  with  a  high  pressure  collisional  plasma  is 
characterized  by  a  decrease  in  conductivity  in  relation  to  the  far-field  plasma  because  the 
anode  acts  as  a  catalyst  to  ion-electron  recombination,  and  hence  is  a  “sink”  for  charged 
particles.  The  charged  particle  density  thus  begins  from  a  small  value  at  the  anode 
surface  and  increases  to  the  value  of  the  far-field  plasma  as  the  distance  away  from  the 
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anode  increases.  This  produces  a  small  scale  region  where  charged  particle  diffusion  is 
driven  by  the  gradients  in  the  plasma  properties.  The  plasma  can  be  in  a  state  of 
ionization  in  this  region  because  of  the  depressed  charged  particle  number  densities.  Due 
to  the  property  gradients  in  this  region,  analogies  can  be  made  to  the  “fluid”  and 
“thermal”  boundary  layers;  hence  the  region  is  called  the  “electric  boundary  layer.”  At 
present,  no  arcjet  model  is  believed  to  capture  this  region  with  sufficient  accuracy  and 
rigor.  Several  near-electrode  models  have  been  developed,  however,  which  attempt  to 
more  completely  capture  the  physics  of  this  region  by  making  simplifications  in 
geometry  and  the  regions  outside  the  electric  boundary  layer  [9-12],  Fig.  1.2  shows  a 
schematic  of  the  two  general  current  attachment  modes  in  a  simplified  planar  anode 
geometry. 


Figure  1.2.  Schematic  of  the  two  attachment  modes  in  a  simplified  planar  geometry 

1.2.  Review  of  Relevant  Studies 

There  have  been  a  large  number  of  studies  performed  on  the  transition  from  the 
diffuse  to  constricted  anode  attachment  modes  for  vacuum  arcs,  but  relatively  few  for  the 
higher  pressure  cases.  This  is  because  of  the  historically  large  interest  in  increasing  the 
reliability  of  triggered  vacuum  gaps  and  vacuum  switches.  It  is  instructive,  however,  to 
examine  these  studies  to  understand  similarities  and  differences  between  the  vacuum  and 
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higher  pressure  cases.  Despite  the  differences  in  operating  pressures,  there  may  be  some 
similarities  in  the  mechanisms  of  mode  transition  between  vacuum  and  high  pressure  arcs 
due  to  phenomena  (such  as  ionization  and  Joule  heating)  common  to  both  cases. 

In  a  vacuum  arc,  the  sources  of  the  interelectrode  plasma  are  the  electrodes 
themselves.  The  term  “vacuum  arc”  is  therefore  a  bit  of  a  misnomer,  since  the  gas  near 
the  electrodes  can  become  collisionally  dominated  [13].  The  current  attachment  process 
is  separated  into  four  different  modes;  the  diffuse  arc,  footpoint,  anode  spot,  and  intense 
arc  modes.  At  low  current  densities,  the  arc  exists  in  the  diffuse  arc  mode.  In  this  mode 
the  anode  is  basically  inert,  and  acts  as  a  collector  of  particles  emitted  from  the  cathode 
(ions,  electrons,  and  atoms).  There  is  no  erosion  at  the  anode  and  the  voltage  remains 
relatively  low  and  quiet.  Raising  the  current  density  may  give  rise  to  the  footpoint  mode, 
which  is  characterized  by  luminous  spots  and  anode  melting.  This  mode  appears  only 
under  certain  conditions,  and  sometimes  may  exist  as  an  interim  mode  between  the 
diffuse  and  anode  spot  modes.  It  is  also  usually  accompanied  by  the  appearance  of  arc 
noise  (voltage  fluctuations)  and  an  increase  in  the  mean  arc  voltage.  These  voltage 
fluctuations  are  associated  with  fluctuations  in  light  intensity  given  off  by  the  arc,  which 
indicate  that  fluctuations  in  the  electron  number  density  and  temperature  [14]  are 
occurring.  The  anode  begins  to  take  an  active  role  in  this  discharge.  The  remaining  two 
modes,  which  are  associated  with  higher  currents,  are  characterized  by  intense  anode 
heating  and  erosion.  The  voltage  is  typically  low  and  quiet  in  these  modes.  The  anode 
acts  as  a  copious  source  of  vapor  and  ions,  and  has  a  temperature  near  the  atmospheric 
boiling  point  of  its  material  [13]. 

Many  mechanisms  have  been  postulated  to  be  responsible  for  the  transition  of  arc 
modes,  as  there  has  not  been  one  mechanism  identified  which  is  exclusively  responsible. 
These  mechanisms  include  magnetic  constriction  in  the  plasma,  as  well  as  fluxes  of 
material  from  the  anode  and  cathode  and  thermal,  geometric,  and  electrical  effects  of  the 
anode.  For  example,  some  researchers  contend  that  a  significant  anode  vapor  pressure, 
which  is  associated  with  higher  anode  temperatures,  is  a  key  requirement  for  transition. 
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Specifically,  Wieckert  and  Egli  [15]  hypothesize  that  plasma  production  from  ionization 
resulting  from  sufficient  neutral  number  densities  near  the  anode  may  provide  the  starting 
point  for  transition  to  an  anode  spot.  Barinov  and  Smirnov  [16]  showed  that  anodes 
composed  of  materials  whose  erosion  and  ionization  characteristics  tended  to  produce 
greater  densities  of  ions  near  the  anode  had  correspondingly  lower  critical  currents  (the 
current  at  which  an  anode  spot  forms)  than  other  electrode  materials.  Others  have  pointed 
out  that  the  increase  in  voltage  noise  and  magnitude  at  the  moment  of  transition  is  due  to 
the  appearance  of  an  ion-  deficient  region  near  the  anode  that  creates  a  significant  voltage 
drop,  which  can  be  on  the  order  of  the  ionization  potential  of  the  plasma,  due  to  the  space 
charge  [13].  The  appearance  of  a  rapidly  oscillating  voltage  at  transition  also  suggests 
the  occurence  of  unstable  phenomena  [17].  Miller  [13]  has  suggested  that  the  most 
probable  explanation  of  anode  spot  formation  is  a  combination  theory,  in  which  many 
effects  contribute  in  some  way. 

While  the  majority  of  research  has  been  of  an  experimental  nature,  there  have  also 
been  a  few  theoretical  studies.  One  notable  attempt  was  that  of  Ecker  [18],  who 
calculated  critical  current  densities  and  anode  temperatures  for  anode-spot  transition.  To 
make  the  calculations  tractable,  Ecker  used  a  very  much  simplified  one-dimensional 
model  which  assumed  collisional  dominance,  local  thermodynamic  equilibrium,  Saha 
equilibrium,  and  neglected  inhomogenities,  conduction  losses,  excitational  energy  losses 
and  the  sheath.  Ecker  found  from  his  model  that  the  transition  was  due  to  an  instability 
which  appears  in  the  arc  near  the  anode,  and  that  a  necessary  but  insufficient  condition 
for  the  instability  is  the  appearance  of  significant  anode  vaporization.  His  results 
exhibited  some  qualitative  agreement  (within  an  order  of  magnitude)  and  some 
differences  with  the  available  experimental  data. 

In  contrast,  there  are  far  fewer  studies  which  discuss  mode  transition  for  higher 
pressure  arcs,  where  a  gas  takes  on  the  role  of  the  interelectrode  plasma.  It  was 
mentioned  in  a  review  paper  by  Kimblin  [19]  that  anode  spot  formation  occurs  at  “short 
electrode  spacings  independent  of  the  arc  current”  in  the  presence  of  an  ambient  gas. 
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This  is  because  the  gas  shields  the  anode  from  the  flow  of  ions  from  the  cathode,  creating 
an  ion-deficient  zone,  which  as  discussed  above  is  said  to  be  a  destabilizing  influence, 
though  many  other  factors  also  have  an  important  effect.  Dinulescu  and  Pfender  [12] 
commented  that  the  addition  of  a  strong  convective  cathode  jet  also  has  a  stabilizing 
affect;  when  electrode  spacing  was  increased,  the  effect  of  the  jet  on  the  anode  was 
reduced,  which  resulted  in  a  mode  transition  at  modest  current  densities. 

Theoretical  studies  on  the  subject  of  atmospheric  arcs  have  been  in  the  form  of 
steady  state  calculations  of  plasma  properties  [9-12].  From  these  studies,  a  better 
understanding  of  the  structure  of  the  anode  sheath  and  non-equilibrium  boundary  layer 
was  obtained.  In  particular,  the  past  work  of  Self  and  Eskin  [9],  Cappelli  [10],  and 
Meeks  and  Cappelli  [11]  have  helped  clarify  the  physical  processes  and  conditions 
distinguishing  the  diffuse  and  constricted  modes  of  current  attachment.  In  those  studies, 
it  was  shown  that  for  an  anode  in  contact  with  a  thermal  plasma,  the  boundary  layer 
voltage  first  increases  with  increasing  current  density,  eventually  reaching  a  local 
maximum.  This  point  defines  the  demarcation  between  a  positive  and  negative  boundary 
layer  impedance,  where  ionization  within  the  electrical  boundary  layer  is  supported  by 
increased  ohmic  heating.  This  current- voltage  layer  was  found  to  be  strongly  influenced 
by  the  external  flowfield.  It  was  speculated  that  the  regions  of  positive  and  negative 
boundary  layer  impedance  correspond  to  the  cases  of  diffuse  and  constricted  modes  of  arc 
attachment.  This  contention,  however,  can  only  be  confirmed  or  denied  by  a  detailed 
unsteady  analysis  of  the  problem  or  by  an  experiment. 

Other  researchers  have  attempted  to  model  the  arc  attachment  process  within  the 
context  of  a  practical  device.  For  example,  Paik  et  al  [20]  have  applied  the  so-called 
Steenbeck  Minimum  Principle  to  a  direct  current  (DC)  discharge  geometry  often 
encountered  in  conventional  arc  plasma  torches.  They  show  that  this  principle,  which  is  a 
special  case  of  minimum  entropy  production,  is  useful  in  determining  the  position  of  the 
anode  arc  root  in  a  plasma  torch.  They  also  find  that  the  axial  distance  from  the  cathode 
tip  to  the  steady  state  arc-root  position  decreases  for  increasing  total  current  due  to  a 
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Okazaki  et  al  [21,22]  have  made  a  serious  attempt  at  theoretically  understanding 
mode  transition  in  an  MHD  generator.  Their  model  identifies  a  parameter,  defined  as  a 
ratio  of  Joule  heat  to  conductive  heat,  which  determines  the  stability  of  the  boundary 
layers  (electrical,  thermal,  and  fluid).  They  conclude  that  a  maximum  value  of  this 
parameter  exists  which,  when  exceeded,  leads  to  instability  and  hence  mode  transition. 
The  analysis  also  relates  stability  characteristics  to  electrode  temperatures,  concluding 
that  lowering  the  anode  temperature,  which  lowers  the  conductivity  of  the  adjacent 
plasma,  reduces  the  critical  current  density  at  which  mode  transition  takes  place.  Their 
model  takes  into  account  ionizational  nonequilibrium  and  inhomogeneity  in  the  plasma 
properties  due  to  the  boundary  layer,  but  also  makes  some  important  simplifications. 
The  plasma,  which  is  formed  from  potassium-seeded  combustion  products,  is  assumed  to 
have  a  perfect  collisional  coupling  between  the  electrons  and  neutrals  (though  L.T.E  is 
not  assumed).  Thus  a  simplified  form  of  the  electron  energy  equation,  which  balances 
inelastic  collisional  energy  losses  to  Joule  heating,  is  used,  ignoring  many  other 
mechanisms  of  electron  energy  transfer  such  as  convection,  conduction,  pressure  work, 
and  energy  dissipation  due  to  ionization. 

It  is  clear  from  this  review  that  the  transition  from  diffuse  to  constricted  current 
modes,  which  is  seen  in  all  regimes  of  plasma  pressures,  can  be  initiated  by  many 
different  mechanisms.  In  the  most  general  case,  many  phenomena,  such  as  three 
dimensionality,  unsteady  behavior,  bulk  fluid  flow,  ion  motion,  ionizational  non¬ 
equilibrium,  electron  energy  conduction,  surface  effects,  etc.  should  be  accounted  for. 
Simplifications  must  be  made  since  this  is  clearly  a  monumental  task.  But  as  the  old 
Taoist  proverb  goes,  “A  journey  of  a  thousand  miles  begins  with  a  single  step.” 

1.3.  Approach 

In  this  thesis  a  theoretical  analysis  of  the  transient  breakdown  phenomenon  in  a 
high  pressure  plasma  is  presented  and  discussed.  The  approach  used  is  to  analytically 
examine  the  stability  characteristics  of  the  electrical  boundary  layer  in  a  simplified  planar 
geometry.  This  particular  study  is  an  extension  of  the  work  of  Self  and  Eskin  [9],  which 
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introduces  a  model  description  for  steady  state  current  transfer  from  a  planar  electrode  to 
a  non-flowing  thermal  plasma.  In  this  previous  work,  steady-state  solutions  of  a 
simplified  set  of  plasma  conservation  equations,  which  describe  the  characteristics  of  the 
one-dimensional  electrical  boundary  layer,  are  obtained.  The  results  of  this  study 
provide  baseline  profiles  which  are  analyzed  for  their  stability  characteristics.  From 
these  results  conditions  which  are  favorable  for  instabilites  to  occur  (and,  we  speculate, 
for  the  mode  of  current  transfer  to  change)  will  be  discussed  and  intuition  about  the 
effects  of  various  physical  mechanisms  on  stability  will  be  developed.  The  model  takes 
into  account  finite-rate  ionization  and  recombination,  the  presence  of  an  ambipolar  field, 
pressure  forces  due  to  gradients  in  number  density  and  temperature,  collisional  energy 
and  momentum  exchange,  pressure  work,  inertia,  and  electron  energy  conduction  and 
convection,  many  of  which  have  been  neglected  in  previous  studies.  Some  important 
effects  such  as  bulk  fluid  motion,  heavy  species  thermal  behavior,  and 
multidimensionality  are  not  included  in  the  model  for  the  sake  of  simplicity,  though  they 
should  be  included  in  future  studies.  Limitations  of  the  model  will  become  apparent  and 
speculation  on  the  effects  of  ignored  processes  will  then  be  given. 
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Chapter  2.  Model  Description 


The  following  1-D  continuum  model  is  based  on  the  equations  derived  by  Eskin 
[9].  These  equations  are  based  on  moments  of  the  Boltzmann  equation  which  describe  the 
behavior  of  electrons  and  ions  in  a  plasma.  As  stated  in  Eskin’ s  thesis,  the  goal  of  this 
model  is  to  examine  the  nature  of  the  plasma-electrode  interactions,  and  not  to  simply 
analyze  one  specific  experiment  or  setup. 

2.1.  Governing  Equations 

The  derivation  of  the  governing  equations  which  we  will  use,  in  their  most  general 
form,  can  be  found  in  Appendix  A  of  Eskin’s  thesis.  It  follows  the  standard  methodology 
of  taking  moments  of  the  Boltzmann  equation  for  the  ions  and  electrons  in  a  weakly 
ionized  plasma  and  making  the  usual  heuristic  descriptions  of  the  collisional  terms.  For 
this  study,  the  inertia  terms  in  the  momentum  equations  were  kept  to  make  the  model 
applicable  to  a  wider  range  of  conditions;  in  addition,  the  unsteady  terms  were  retained 
because  of  temporal  nature  of  the  investigation.  The  dependent  variables  are  the  electron 
and  ion  number  densities,  the  electron  and  ion  velocities,  the  electric  field,  and  the  electron 
temperature. 

The  set  of  equations  describing  the  behavior  of  the  electrons  and  ions  are  as 
follows: 

Ion  (j=i)and  Electron  (j=e)  Conservation 


+  (fijVj  )  =  he=  nen<x(Te )  -  n,nt2p(  Te ) 


(2.1.1) 
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Here,  n;  and  Vj  represent  the  number  density  and  velocity  of  species  j  (either  electrons  or 
ions),  n  represents  the  background  neutral  species  density,  he  is  the  volumetric  rate  of 
ionization  (for  both  the  electrons  and  ions),  and  Te  is  the  electron  temperature.  We 
assume  that  electron  collision-induced  ionization  and  recombination  dominate  in  the 
ionization  and  recombination  process,  and  that  P(Te)is  the  three  body  recombination 
coefficient  given  by  Hinnov  and  Hirschberg  [23]: 

p(Te )  =  1.09  x  1 O"20 T,~ 2  y  (2. 1 .2) 


For  a  weakly-ionized  plasma,  the  ionization  rate  coefficient,  cx(Te),  is  related  to  P(Te),  the 
neutral  temperature,  T,  and  the  pressure,  P,  through  detailed  balance: 


(2.1.3) 


Ion  (j=i)  and  Electron  (j=e)  Momentum  Conservation 


dV: 


n]m]  -y-  +  n  -m,  {v ■  ■  VjV;  =  —VP,  +  z}en}E  -  m-np - 


(2.1.4) 


where  ze  =  -1,  z,  =  +1 ,  e  is  the  magnitude  of  the  electron  charge,  Pe  =  nekTe  and 
Pi  =  n;kT  are  the  electron  and  ion  pressures  respectively,  and  v;  is  the  electron  or  ion 
collision  frequency  with  the  background  gas  (related  to  the  corresponding  electron  and  ion 
mobility  /ij  =  e/mjVj).  We  have  neglected  collisions  between  charged  particles  because  of 
the  assumed  weakly  ionized  condition,  the  momentum  transferred  to  the  bath  of  electrons 
and  ions  due  to  the  finite  rate  of  ionization/recombination,  and  shear  stresses. 
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Electron  Energy 

3  dT  3  -  _  _  m 

2  n^ir + 2  V7; + v  •  ^ + ■  m*v”r‘  • = 3<5*  v”n*k(T  ~T*y £ A 

(2.1.5) 

Here  m  is  the  mass  of  the  background  neutral  species,  &  is  the  nonelastic  electron  energy 
loss  factor,  which  is  much  greater  than  unity  only  for  molecular  plasmas  [24],  and  £j  is  the 
ionization  energy  of  the  neutral  species.  We  have  neglected  work  on  electrons  by  shear 
stresses. 

Poisson’s  Equation 


V  i  =  —  (n{  -ne) 

£ 

co 


(2.1.6) 


where  eo  is  the  permitivity  of  free  space  ( 8.854  x  1 0  12  Farads/m) . 

The  heavy  species  (ions,  neutrals)  are  assumed  to  be  efficiently  coupled,  and 
therefore  share  the  same  uniform  temperature,  T,  which  removes  the  need  for  a  global 
energy  equation.  The  total  pressure,  P,  is  determined  from  the  equation  of  state: 

P  =  (ne  +  n)kT  +  nekTe  ~  nkT  (2. 1 .7) 

The  electron  and  ion  momentum  equations  are  used  to  calculate  the  electron  and  ion 
velocities.  The  conduction  current  density  can  be  calculated  through  the  ancillary  relation: 


J  =  Je  +  Ji=ene(vi-Ve) 


(2.1.8) 
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2.2.  Assumptions 

The  additional  assumptions  which  we  use  to  further  simplify  the  model  to  make  the 
calculations  more  tractable  are  given  as  follows.  The  independent  spatial  variable,  x, 
describes  the  distance  from  the  anode. 

1 .  A  one  dimensional  system  with  a  plane  electrode  contacting  a  semi¬ 
infinite  plasma  in  the  half  space  0  <  x  <  «>,  as  shown  in  Fig  2. 1 . 

Far-field  Plasma 


Current 


Current 


Planar  Anode 


Figure  2.1 .  Schematic  of  the  simplified,  1-D  anode  model. 

2.  Ions  and  background  neutral  particles  are  at  the  same  uniform  temperature. 

3.  There  is  no  external  magnetic  field,  and  such  fields  due  to  plasma  currents 
are  neglected. 

4.  A  net  current  to  or  from  the  electrode  is  allowed,  and  is  assumed  to  be  due 
to  an  applied  potential  between  the  electrode  and  the  distant  plasma. 

5.  The  plasma  is  weakly  ionized  and  quasineutral. 

6.  There  is  no  bulk  fluid  flow. 
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By  invoking  quasineutrality,  it  is  implied  that  the  scale  length  of  the  phenomena  of 
interest  is  much  larger  than  the  Debye  length,  and  that  the  frequency  of  any  unsteady 
behavior  is  smaller  than  the  plasma  frequency.  The  requirement  that  heavy  particles  are  at 
a  constant  uniform  temperature  is  akin  to  assuming  that  the  specific  heat  or  conduction 
coefficient  of  the  heavy  species  is  infinitely  large. 

2.3.  Non-Dimensionalization 

To  keep  the  model  as  general  and  simple  as  possible,  all  variables  are  non- 
dimensionalized  by  dividing  them  by  reference  parameters  which  correspond  to  the  distant 
plasma  conditions  for  the  initially  steady  state  one-dimensional  electrical  boundary  layer. 
The  reference  parameters  are  defined  in  Table  2. 1 . 

The  variables  are  non-dimensionalized  as  follows: 


(2.3.1) 


U,  = 


U; 


Ir  / 


where  U  is  the  x-component  of  the  velocity  vector  V  . 


-  ri 

ij  =  — - ,  where  i  is  the  x-component  of  the  flux  vector  JT . 
Fr 


The  non-dimensional  set  of  equations  which  result  from  the  substitution  of  these  non- 
dimensional  variables  become: 
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Table  2.1.  Reference  Parameters  for  Non-dimensionalization 


nD 


n„  =  P(T„  H 


1  2  Da 

4=. 

mTJn » 

2D„nR 

rR  = 

— -  mUR 

kT 

Er  ~ 

eR 

elR 

Jr  = 

GrEr 

’’-pfcH 

where 

^R=nRe(fJi+fie) 

D  _  A  a,  +  Aa 
“  /4  +  A 


D 


Z>  = 


kTCfJle 

e 

kTe  fa 


:  farfield,  zero  current,  equilibrium  electron  number  density 
:  the  zero  current,  distant  plasma  electron  temperature  (=7) 
:  the  recombination  rate  in  the  distant  plasma 

:  the  electron-ion  recombination  scale  length 

:  the  characteristic  ambipolar  flux 

:  the  characteristic  ambipolar  field 
:  the  characteristic  current  density 
:  the  distant  plasma  characteristic  recombination  time 

:  the  distant  plasma  electrical  conductivity 
:  the  distant  plasma  ambipolar  diffusion  coefficient 

:  the  distant  plasma  electron  free-diffusion  coefficient 
:  the  distant  plasma  ion  free-diffusion  coefficient 


Ion  and  Electron  Conservation 


dt  dx 


Te2  exp 


1 


(2.3.2) 
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where,  for  convenience,  we  have  substituted  m  =  ne  from  quasineutrality. 

Subtracting  the  electron  conservation  equation  from  the  ion  conservation  equation 
gives  a  simple  relationship  between  the  ion  and  electron  flux: 

|(r,-r.)=o  (2.3.3) 


or  J  -  Jfarfieid  -  constant.  This  is  a  limitation  of  the  1-D  quasineutral  model  which  does 
not  allow  us  to  investigate  the  more  general  case  of  a  non-zero  perturbation  in  current 
density.  It  can  be  removed  by  expanding  the  model  to  include  two-dimensional  effects. 

Electron  Momentum  Conservation 


m. 


(j  +  dKC 


m 


'dU. 


dt 


■r  +  U . 


dJL 

dx 


=  -n,E  - 


-  dneTe  4ji 


dx  1  +  ji 


(2.3.4) 


Ion  Momentum  Conservation 


fdU;  -  dU; 


(l  +  v)vjR 


n. 


dt 


d-  +  U; 


dx 


=  n„E  — 


dn„ 


dx  1  +  fi 


4  E 


(2.3.5) 


Electron  Energy  Conservation 


3  _  dT  3  _  —  dT  _  —  dU  3  1  +  fl  d 
2  e  di  2  e  e  dx  e  e  dx  5  ii  dx 


( 


-rdC 
neTe  — 


dx 


4fi 

1  +  p 


reUe=36ne{l-Te)-£,h 


(2.3.6) 
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The  dimensionless  parameters  which  we  introduce  in  the  electron  energy  equation 
ar &9  ■=  8HtRe  / mifle,  the  collision  coupling  factor,  ji  =  fii  / pte,  the  ratio  of  ion  and 
electron  mobility,  and  'el-eI/  kTe  ,  the  non-dimensional  ionization  energy 

2.4.  Stability  Analysis 

The  stability  of  the  electrical  boundary  layer  is  analyzed  by  modeling  the  time 
dependent  effect  of  small  oscillating  perturbations  in  the  plasma  properties  ne,  Ue,  Ut,  E,  Te 
which  are  superimposed  on  a  set  of  steady  state  initial  conditions.  In  keeping  with  our  1- 
D  geometry,  we  assume  that  the  perturbations  are  planar  and  hence  propagate  and  vary  in 
the  x-direction  only,  perpendicular  to  the  planar  anode.  Mathematically,  each  variable  is 
decomposed  into  a  component  representing  the  steady  solution  and  one  representing  the 
small  perturbation  as  follows, 

f{x,t)  =  fjx)  +  f'(x,t)  (2.4.1) 

where  “f  denotes  a  plasma  property  such  as  ne.  Ue,  Ut,  E,  Te,  and 

f'{x,t)  =  f{x)e-m  (2.4.2) 

corresponds  to  a  spatially  varying  temporally  damped  or  growing  oscillating  perturbation. 
Here  co,  which  is  in  general  complex,  represents  both  the  frequency  of  oscillation  (real 
component)  and  growth  (positive  imaginary  component)  or  decay  (negative  imaginary 

component)  rate.  The  eigenfunction  /( x) ,  which  in  general  is  also  complex,  describes  the 
spatial  amplitude  as  well  as  the  phase  relationships  between  the  set  of  perturbation 
variables.  There  is  no  loss  in  generality  by  the  use  of  this  form  for  the  disturbances 
because  of  the  fact  that  arbirtrary  perturbations  to  linear  systems,  into  which  we  can 
decompose  our  problem  with  a  small  amplitude  approximation,  can  be  decomposed  into  a 
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sum  of  fourier  components.  Each  fourier  component  is  represented  by  its  frequency  a  and 
can  be  analyzed  independently  from  the  others. 

The  focus  of  this  study  is  to  find  situations  where  the  imaginary  component  of  ox 
OX,  is  positive.  This  denotes  a  perturbation  which  is  growing  exponentially  in  time,  and 
corresponds  to  a  linear  instability.  Linear  stability  theory  can  determine  whether  there  is 
an  initial  growth  of  a  small  disturbance,  but  cannot  determine  the  behavior  of  that 
disturbance  once  its  amplitude  becomes  large  enough  for  non-linear  effects  to  become 
important.  As  a  result,  linear  stability  theory  does  not  guarantee  that  an  unbounded 
growth  is  going  to  actually  occur;  however,  because  of  the  simplifications  provided  by 
using  linear  theory,  it  is  often  employed  as  a  first  step  for  gaining  insight  on  the  stability 
characteristics  of  physical  situations. 

To  obtain  the  stability  equations,  equation  2.4.1  is  substituted  into  the  general 
equations  2.3.2  -  2.3.6  for  each  plasma  property  ( ne ,  Ue,  Uu  E,  Te).  Because  the  analysis 
is  applied  to  known  steady  state  solutions,  terms  which  also  appear  in  the  steady  state 
equations  drop  out.  With  the  assumption  that  the  perturbed  quantities  are  very  small 
compared  to  the  steady  state  quantities,  products  of  perturbed  quantities  can  also  be 
eliminated;  the  set  is  thus  linearized.  The  remaining  terms  then  form  the  following  set  of 
linear  ordinary  differential  equations: 

Particle  Conservation 


^  dr  -  - 

-io)ne  +  — f-  =  Te~J  exp 
dx 
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(2.4.3) 


~  dT;  ± 
■icme+  —  =  ne 
dx 


(2.4.4) 
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Ion  Momentimi  Conservation 


(l  +  p)vjR 


^  _  dll,  „  dU,  ^  _  dU; 

-icon.  U.  +U:  —^n.+n.  —^rU;+K  Ut 


^  dne  4 

=  n.  E  +  Es,ne  -  ——  — - 1  n 

"  e  dx  1  +  pi 


dx  "*  dx  dx 

&,+v,i.) 


(2.4.5) 


Electron  Momentum  Conservation 
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(2.4.6) 


Electron  Energy  Conservation 
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(2.4.7) 


This  set  of  equations  describes  the  spatial  variation  of  the  perturbations  f{x  )  with 

distance  from  the  electrode,  and  are  subject  to  appropriate  homogeneous  boundary 
conditions  at  the  far-field  plasma  and  anode  surface.  These  boundary  conditions  will  be 
discussed  in  the  next  section.  The  actual  solution  to  this  set  of  equations  is  obtained  by 
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solving  a  generalized  eigenvalue  problem,  where  co  is  the  eigenvalue  and  the  perturbation 
amplitudes  We,Ue,UiyE ,Te  form  the  eigenvector. 

2.5.  Boundary  Conditions 

Equations  2.3.2  -  2.3.6,  without  their  unsteady  terms,  are  a  sixth-order,  coupled, 
nonlinear  set  of  partial  differential  equations.  Six  boundary  conditions  are  therefore 
required  for  their  steady-state  solution,  with  5  initial  conditions  additionally  required  for 
the  unsteady  problem.  For  the  stability  analysis,  the  steady  state  solution  acts  as  the  initial 
condition.  Because  of  this,  the  solution  of  the  steady  problem  is  considered  first. 

In  the  farfield  plasma,  which  we  assume  to  be  uniform,  the  four  boundary 
conditions  are  found  via  the  solution  of  the  following  set  of  algebraic  equations: 


(  ,  w 
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exp 
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\  Ie~  J) 

—  1  +  /t  E  J 

T  =1  + - — — — 

4fi  39ne_ 


l+E 
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r. 


i+fi 

4/1 


J 


(2.5.1) 

(2.5.2) 

(2.5.3) 

(2.5.4) 

(2.5.5) 


The  assumption  that  Saha  equilibrium  exists  in  the  far-field  region  is  the  most 
convenient  one  possible,  because  it  allows  the  far-field  region  to  be  uniform.  In  some 
devices,  such  as  the  arcjet  thruster,  the  region  far  from  the  anode  is  not  in  ionizational 
equilibrium  and  hence  is  not  uniform.  For  these  cases,  additional  application  specific 
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information  is  needed  to  formulate  the  far-field  boundary  condition.  However,  in  order  to 
keep  the  model  as  general  and  simple  as  possible,  the  assumption  that  equilibrium  and 
uniform  conditions  exist  in  the  farfield  is  utilized  in  this  study. 

For  the  stability  analysis,  all  perturbations  are  set  to  zero  in  the  far-field.  This 
enforces  the  condition  that  the  disturbances  are  not  felt  far  away  from  the  anode.  In  other 
words,  the  disturbances  are  assumed  to  originate  in  the  electric  boundary  layer  and  die  off 
away  from  the  anode.  Another  possibilty  is  to  assume  a  constant  amplitude  oscillation  in 
the  far-field  region;  this,  however,  can  imply  that  disturbances  are  originating  away  from 
the  anode  and  propagating  towards  it.  There  is  also  the  possibility  that  these  two 
boundary  conditions  represent  the  same  disturbance.  For  example,  a  perturbation  whose 
amplitude  settles  to  a  zero  gradient  in  the  far-field  region  may  also  be  settling  towards  a 
zero  value.  For  this  study,  it  has  been  observed  that  both  boundary  conditions  produce 
the  same  results. 

2.6.  Anode  Surface  Boundary  Condition 

Two  additional  boundary  conditions  are  required  at  the  anode  for  the  steady  state 
calculations.  There  have  been  many  methods  devised  to  model  the  way  the  plasma 
interacts  with  the  anode.  Some  of  these  have  been  employed  because  of  their 
mathematical  simplicity.  Examples  of  this  are  specifying  the  electron  and  ion  number 
densities  or  fluxes  at  the  anode.  Although  these  representations  may  allow  mathematical 
solutions  for  the  problem,  care  must  be  exercised  in  the  interpretation  of  results.  Recently 
there  have  been  attempts  to  use  boundary  conditions  which  possess  more  basis  in  physical 
reasoning.  Examples  of  these  methods  are  employing  a  unity  sticking  coefficient  for 
electrons,  modeling  the  emission  characteristics  of  the  electrode,  and  kinetic  approaches. 
A  more  exhaustive  review  of  the  various  methods  which  researchers  utilized  can  be  found 
in  Meeks’  thesis  [11]. 
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The  difficulty  in  modeling  the  surface  boundary  condition  is  that  the  anode 
interacts  in  many  different  ways  with  the  plasma.  Meeks  has  come  up  with  a  fairly  broad 
list  of  these  interactions  in  her  thesis,  which  include  thermionic  and  secondary  electron 
emission,  ion-surface  charge  exchange,  neutral-surface  charge  exchange,  electron 
adsorption  and  the  ability  of  the  elctrode  to  capture  or  release  “bulk”  electrons 
contributing  to  current  flow  (where  “bulk”  refers  to  the  solid  electrode  material). 
However,  because  the  individual  reaction  rates  for  all  these  mechanisms  are  unknown,  a 
global  mechanism  which  takes  the  place  of  all  these  separate  effects  is  used.  This  is 
implemented  with  an  expression  for  the  global  surface  reaction  rate,  which  is  postulated  to 
have  a  first-order  dependence  on  both  the  electron  and  ion  densities  as  follows, 

T;(x  =  0)  =  -ksin2e{x  =  0)  (2.6.1) 

where  the  quasineutral  approximation  has  been  applied.  kSi  is  the  reaction  rate  for  the 
surface  reaction  represented  by  : 

A+  +  e'=>A  (2.6.2) 

Although  ksi  is  still  unknown,  it  is  only  a  single  parameter  and  is  thus  much  less 
cumbersome  to  deal  with. 

For  the  steady  state  solutions,  the  condition  that 

— f-  (x  =  0)  =  0  (2.6.3) 

dx 

is  used.  This  is  the  default  boundary  condition  used  in  some  arcjet  models  [4],  and  is  also 
used  in  this  study  for  consistency. 
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2.7.  Approach  to  Solutions 

The  basic  approach  taken  in  this  study  is  to: 

1.  Establish  a  baseline  steady  state  solution  methodology  and  results.  The  characteristics 
of  the  electrical  boundary  layer  will  be  presented  and  reviewed. 

2.  Perform  a  simplified  stability  analysis  based  on  a  hypothetical  uniform  underdense 
plasma.  The  problem  is  greatly  simplified  while  providing  initial  insight  into  the 
problem. 

3.  Extend  the  stability  analysis  to  take  into  account  the  non-uniformity  of  the  steady  state 
electrical  boundary  layer.  Here  the  most  general  case  of  this  study  is  presented. 
Insights  from  the  previous  two  sections  are  used  and  new  information  and  insights  are 
discussed. 


Chapter  3. 

Steady  Behavior  of  the  Electric  Boundary  Layer 
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We  begin  the  study  with  the  solution  of  the  steady  state  equations  which  describe 
the  background  plasma  in  the  stability  analysis.  The  steady  state  equations  are  equations 
2.3.2  -  2.3.6  without,  of  course,  the  unsteady  terms. 

3.1.  Non-dimensional  coefficients 


Some  insight  into  the  character  of  these  equations  can  be  obtained  by  examining 
the  constant  coefficients  in  the  equations.  For  example,  by  comparing  the  coefficient  of 
the  inertia  term  in  the  ion  momentum  equation  with  the  coefficient  of  the  collisional 
term,  we  see  that  if  the  product  of  tr,  the  characteristic  recombination  time,  and  vin,  the 
ion-neutral  collision  frequency,  is  much  larger  than  order  1,  then  the  inertia  term  is 
insignificant  in  comparison  with  the  collisional  exchange  term.  Hence,  under  the  higher 
pressure  conditions  for  which  this  is  true,  collisional  forces  dominate.  A  similar 
argument  can  be  made  for  the  electron  momentum  equation,  where  we  have  the 
following  condition  for  collisional  dominance: 


me  l 
m,  vm 


«jl 


(3.1.1) 


The  parameter  6  determines  how  “collisionally  coupled”  the  electrons  are  with  the 
heavy  species.  Large  values  of  6  denote  situations  of  approximate  thermal  equilibrium 
between  the  electrons  and  heavy  species,  whereas  small  values  denote  less  energy  exchange 
and  hence  greater  differences  between  the  electron  and  heavy  species’  temperatures. 
Additionally,  because  of  the  assumption  of  a  uniform  heavy  particle  temperature,  very  high 
values  of  6  result  in  an  isothermal  condition  for  the  electrons  in  this  model. 
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3.2.  Baseline  conditions 

To  get  an  idea  of  what  the  magnitudes  of  the  parameters  9,  p  are  in  a  real 

plasma  of  interest,  we  begin  by  examining  the  steady  conditions  which  will  become  the 
basis  for  our  stability  analyses  in  this  study.  The  situation  we  choose  to  model  is  the 
current  attachment  zone  in  a  typical  1-kw  arcjet  thruster  running  on  hydrogen.  We  again 
note  that  only  a  full  multi-dimensional  model  can  accurately  capture  the  arcjet  boundary 
layer  profile,  and  that  the  results  given  here  are  intended  to  give  qualitative  insight  and 
approximate  behavior.  Specifically,  we  additionally  simplify  the  model  by  assuming  that 
the  plasma  is  10%  dissociated  (i.e.  the  partial  pressure  of  the  dissociated  atoms  is  10%  of 
the  total  pressure),  and  that  the  primary  collisions  are  between  electrons  and  hydrogen 
(HL,)  molecules  and  dissociated  hydrogen  ions  (H+)  and  hydrogen  molecules  only. 
Charged  particle  collisions  are  thus  ignored.  Collisional  cross  section  values  are  roughly 
estimated  from  the  representitive  data  presented  in  Mitchner  and  Kruger  [24]. 
Collisional  frequency  calculations  involving  the  electrons  are  also  based  on  the  bulk 
temperature,  which  is  a  known  parameter,  rather  than  electron  temperature,  which  must 
be  calculated.  The  values  for  the  parameters  used  are  given  in  Table  3.1. 

Equations  2.3.2  -  2.3.6  can  now  be  solved  given  these  parameters.  The  solution 
method  we  utilize  in  this  study  is  based  on  a  simple  finite  difference  scheme  which  uses 
one-sided  differences  for  first  order  derivatives  and  central  differences  for  the  second 
order  derivatives.  The  discretized  matrix  is  solved  using  a  modified  Newton-Raphson 
method  present  in  the  subroutine  TWOPNT  [26].  The  results  of  the  calculation  provide 
the  spatial  profiles  of  ne ,  Ue,  U;,  E  ,  Te,  and  h  for  a  given  current  J .  For  this  study, 

results  are  given  for  three  non-dimensional  current  densities  of  0.4,  40,  and  400.  These 
correspond  approximately  to  dimensional  current  densities  of  10\  104,  and  10s  A/m2, 
respectively. 

We  do  not  include  the  J=  4  case  because  the  boundary  layer  displays  some  unphysical 
characteristics  in  the  current  range  from  7  =  2  to  7  =  20  for  these  conditions.  The 
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Table  3.1.  Property  Values  and  Conditions  for  Baseline  “Arcjet”  Case 


Gas 

Hydrogen 

Ionization  energy 

13.6  ev 

Pressure,  P 

.5  atm 

Bulk  Temperature,  T 

6000  K 

Qen 

1019  m2 

Qin 

O 

00 

3 

2.94xl010  s’1 

8.41xl09  s'1 

rrif 

9.1095xl0'31  kg 

mi 

1.6744x1 O’27  kg  =  m} 

K 

10IS  cm4  mole'1  s1 

4 

5 

K 

1.33xl0'6  m 

ner 

1.61xl019  m'3 

4 

.0356  s 

4 

.0289  m 

er 

17.9  N/C 

hR 

4.52xl020  m'3sl 

e 

2.85xl06 

£i 

26.3 

ft 

.0019 

range  of  current  densities  for  which  this  condition  occurs  depends  upon  the  value  of  ksi. 
At  the  upper  and  lower  limits  of  the  range  the  ion  flux  is  equal  to  zero  at  the  anode 
boundary.  At  current  densities  within  the  saturation  range,  the  calculated  solution 
requires  that  the  ion  flux  at  the  anode  be  positive.  This  phenomena,  identified  as  current 
saturation,  has  been  discussed  in  previous  studies  of  steady  state  electric  boundary  layer 
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behavior  [10,11],  It  occurs  for  current  densities  in  which  the  charged  particles  are  not 
produced  in  sufficient  quantities  to  satisfy  a  negative  flux  of  ions  to  the  anode.  Because 
anodes  are  very  poor  ion  emitters,  a  positive  ion  flux  represents  an  unrealistic  situation. 
The  range  of  saturation  is  also  affected  by  ep  the  ionization  energy  of  the  gas.  Gases 
with  lower  ionization  energies  tend  to  more  easily  supply  the  required  ions  to  the  anode 
and  thus  have  smaller  ranges  of  current  saturation  than  those  with  higher  ionization 
energies.  We  should  also  note  that  it  is  possible  in  such  conditions  for  this  phenomena  to 
not  appear  at  all. 

3.3.  Steady  state  profiles 

The  steady  state  profiles  for  the  three  current  density  cases  are  presented  in  Figs. 
3.1  -  3.8.  To  maintain  consistency  between  these  figures,  the  1  =  400  ,  40,  and  0.4 
profiles  are  presented  as  a  dashed  line,  a  solid  line,  and  a  dotted  line  respectively.  In 
order  to  show  the  profiles  on  the  same  plot,  some  figures,  namely  the  electron  number 
density  and  flux  profiles,  show  profiles  which  have  been  normalized  to  their  far-field 
values.  Some  curves  are  also  displayed  on  specific  length  scales  to  better  examine  the 
particular  regions  of  interest. 

We  begin  the  discussion  with  the  results  of  the  7=  0.4  case.  From  Fig.  3.2, 
which  shows  the  electron  temperature  profile,  we  see  that  the  electrons  are  very  nearly 
equilibrated  with  the  heavy  species.  In  fact,  the  electron  temperature  is  actually  slightly 
below  the  heavy  species  temperature  near  the  anode.  Fig.  3.3,  which  shows  the  electric 
field  profiles,  reveals  that  the  field  is  actually  negative  in  the  boundary  layer  and  hence 
impeding  the  flux  of  electrons  toward  the  anode.  Because  of  this,  a  “Joule  cooling”  of 
the  electrons  is  occuring  near  the  anode.  At  this  low  current  density,  the  ambipolar  field 
must  impede  the  flow  of  electrons  and  ions  because  of  the  strong  driving  force  provided 
by  the  pressure  gradient.  But  despite  these  observations,  the  boundary  layer  remains  in  a 
state  of  ionization  (though  the  ionization  rates  are  very  low  compared  to  the  other  two 


Steady  Behavior  of  the  Electric  Boundary  Layer 


27 


Figure  3.1.  Profiles  of  the  non-dimensional  electron  number  density  for  J  =  400  (dashed),  40  (solid), 
and  0.4  (dotted).  The  curves  are  normalized  by  the  farfield  values,  which  are  4.51,  1.24,  and  1.00 
respectively. 


Figure  3.2.  Profiles  of  the  non-dimensional  electron  temperature  for  J  =  400  (dashed),  40  (solid), 
and  0.4  (dotted). 
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current  density  cases)  as  shown  by  Figs.  3.4  and  3.5,  and  as  required  by  the  three  body 
recombination  occurring  at  the  anode. 


For  higher  current  densities,  such  as  J  =  40  and  400,  we  see  that  the  boundary 
layer  exhibits  additional  characteristics.  Most  importantly,  there  exists  a  much  stronger 
thermal  non-equilibrium  between  the  electrons  and  heavy  particles.  It  is  clear  that  for 
these  particular  cases  there  is  significant  heating  of  the  electrons  near  the  anode.  This  is 
due  to  the  reduced  conductivity  caused  by  the  loss  of  electrons  and  ions  at  the  anode 
wall,  which  in  turn  forces  the  electric  field  to  rise  in  order  to  maintain  the  required 
current  density.  The  elevated  electron  temperatures  contribute  heavily  to  the  increase  in 
charged  particle  production  for  these  higher  current  cases.  As  shown  in  Fig.  3.5,  the 
region  of  ionization  moves  very  close  to  the  anode  wall  because  of  the  combination  of 
maximum  electron  temperatures  and  minimum  charged  particle  number  densities.  This 
is  in  contrast  to  the  7  =  0.4  case,  where  ionization  is  due  to  the  depressed  charged 
particle  number  densities  alone. 

3.4.  The  Assumption  of  Quasineutrality 

As  discussed  in  the  previous  chapter,  the  assumption  of  quasineutrality  removes 
the  ability  to  resolve  the  sheath  in  this  model.  Let  us  take  a  moment  to  comment  on  this 
issue.  The  equation  we  have  left  out  is  Poisson’s  equation: 


dE  _  n;  -  ne 
dx  £0 


(3.4.1) 


which,  non-dimensionalized,  appears  as 
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where  £=Xd/Ir,  the  ratio  of  the  debye  length  XD  =  ^e0kTff  /n^  to  the  recombination 

scale  length.  When  this  equation  is  ignored,  it  is  assumed  that  the  electric  field  gradient 
is  not  great  enough  to  warrant  a  significant  difference  between  nt  and  nc.  Another  way  to 
write  this  is 

2  dE  _ 

e2  —  «ne  (3.4.3) 

ax 

If  this  inequality  is  true,  then  the  difference  in  nt  and  ne  necessary  to  maintain  the  electric 
field  gradient  is  very  small  compared  to  the  actual  value  of  ne.  For  the  conditions 
described  in  Table  3.1,  we  find  that  £  =  4.6xl0'5.  Plugging  simple  order  of  magnitude 
estimates  for  the  field  gradient  shows  that  this  condition  is  satisfied  to  a  reasonable 
degree.  Therefore,  under  these  conditions,  the  assumption  of  quasineutrality  is  a 
reasonable  simplification. 
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We  begin  the  stability  analysis  by  examining  the  simple  situation  of  an  initially 
uniform  background  plasma.  With  this  simplification  much  can  be  learned  about  the 
effects  of  the  physical  mechanisms  which  are  believed  to  cause  or  hinder  instability 
without  dealing  with  the  additional  complications  associated  with  a  spatially  varying 
plasma. 

4.1.  Mathematical  and  physical  preliminaries 

The  assumption  of  an  initially  uniform,  unbounded  plasma  reduces  the 
mathematical  model  to  a  homogeneous  set  of  linear  algebraic  equations,  rather  than  a  set 
of  homogeneous  ordinary  differential  equations.  Because  of  this  uniformity,  an 
additional  simplification  can  be  made  with  regard  to  the  form  of  the  perturbations 
described  by  eqn.  2.4.2.  For  this  special  case,  we  may  assume  that 

Kx)  =  fe*  (4.1.1) 

where  f  is  a  constant  and  k  is  real-valued  and  denotes  the  non-dimensional  wave 
number.  The  perturbation  thus  exhibits  a  pure  sinusoidally  varying  spatial  character. 
When  combined  with  the  oscillatory  temporal  behavior,  the  perturbation  is  interpreted  as 
a  traveling  plane  wave  which  can  grow  or  decay  in  time.  Even  with  these 
simplifications,  the  solution  is  still  general  because  of  the  fact  that  an  arbitrary  1-D 
perturbation  can  be  decomposed  into  a  sum  of  fourier  components  which  have  the 
traveling  waveform  of 
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To  analyze  the  problem  for  stability,  we  solve  for  the  dispersion  relationships  of 
the  waves,  which  involves  solving  a  generalized  eigenvalue  problem  to  find  the 
relationship  between  co  and  k  .  For  a  specified  wavenumber  k  ,  we  specifically  search 
for  conditions  where  the  imaginary  component  of  co  ,  ,  which  is  interpreted  as  the 

growth  rate,  is  positive,  indicating  a  temporal  instability.  The  focus  of  this  study  will 
only  be  on  temporal  instabilities,  as  opposed  to  spatial  instabilities,  which  involve  the 
determination  of  the  spatial  growth  rates  Im(  k )  for  a  specified  real  frequency  co  . 

We  suspect  that  the  spatially  non-uniform  electrical  boundary  layer,  which  was 
shown  in  the  previous  chapter  to  be  in  a  state  of  ionizational  nonequilibrium,  is  more 
susceptible  to  instability  than  the  far-field  equilibrium  plasma.  To  investigate  the  effect 
of  ionizational  non-equilibrium,  we  allow  the  uniform  plasma  to  be  underdense,  meaning 
at  an  electron  number  density  lower  than  that  dictated  by  local  thermodynamic 
equilibrium  at  the  local  electron  number  density  and  pressure.  We  therefore  impose  the 
existence  of  a  constant  ion/electron  sink  on  the  problem  which,  although  unphysical, 
allows  us  to  retain  the  simplicity  of  the  uniform  analysis  while  including  many  important 
physical  mechanisms. 

4.2.  Steady  behavior 

We  begin  the  analysis  with  the  solution  of  the  steady  state  equations  which 
describe  the  background  plasma.  Because  of  the  uniformity  of  the  plasma,  the  equations 
reduce  to  simple  algebraic  relationships  for  ion  velocity,  electron  velocity,  electric  field, 
and  electron  temperature.  Although  we  ignore  the  steady  ion  and  electron  continuity 
equations  and  hypothetically  model  the  presence  of  a  charged  particle  sink  and  neutral 
particle  source,  we  do  take  into  account  the  production  of  charged  particles  due  to 
perturbations  in  electron  number  density  and  temperature  by  including  charged  particle 
conservation  in  the  set  of  disturbance  equations. 
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Without  the  steady  state  continuity  equations,  we  are  left  with  the  electron 
number  density,  in  addition  to  the  current  density,  as  a  free  parameter.  We  are  therefore 
free  to  choose  the  value  of  the  electron  number  density  and  hence  artificially  set  the 
initial  steady  state  rate  of  ionization.  The  steady  state,  uniform  equations  used  are  merely 
eqns  2.3.3  -  2.3.6  (minus  continuity)  without  the  unsteady  and  gradient  terms: 


J(  l+/t) 

(4.2.3) 

4  nefi 

7(i +a0 

(4.2.4) 
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/  ^  J 

(4.2.6) 


Figs.  4.1  and  4.2  give  the  variation  in  these  computed  steady  state  plasma  properties  as  a 
function  of  the  non-dimensional  electron  number  density  for  a  specified  non-dimensional 
current  of  /  =  363  (/=  lxlO5  A/m2).  The  remaining  properties  are  identical  to  those 
listed  in  Table  3.1. 

It  is  apparent  from  Figs.  4.1  and  4.2  that  as  the  electron  number  density  (and 
hence  conductivity)  decreases,  the  electron  temperature,  electron  and  ion  speed,  and 
electric  field  increase. 


4.3.  Perturbation  equations 

With  the  simplifications  provided  by  the  assumption  of  an  initially  uniform 
plasma,  the  general  perturbation  equations,  eqns.  2.4.3  -  2.4.7,  are  reduced  to  the 
following  set  of  linear  algebraic  equations: 
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Figure  4.1 .  The  steady,  non- 
dimensional  electron  (dashed)  and  ion 
(solid)  velocities  vs.  the  electron  number 
density.  The  curves  are  normalized  by 
10s  and  100,  respectively. 


Figure  4.2.  The  steady,  non- 
dimensional  electric  field  (solid, 
normalized  by  250)  and  electron 
temperature  (dashed)  vs.  the  electron 
number  density. 
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These  equations  represent  a  balance  of  stabilizing  (damping)  terms,  destabilizing 
(growth)  terms,  and  terms  whose  effect  on  stability  is  not  obvious  by  inspection.  The 
first  two  types  of  terms  are  those  for  which  the  perturbed  property  (electron  number 
density,  temperature,  etc)  is  the  same  as  that  of  the  unsteady  term  (the  term  which 
contains  of.  For  example,  consider  the  following  linear  ordinary  differential  equation: 

dy 

dt  ay  (4-3.6) 


The  solution  to  this  o.d.e  has  the  form 

y  =  e±of  +  c  (4.3.7) 

which  denotes  either  an  exponentially  decaying  (-)  or  growing  (+)  temporal  solution.  If 
we  substitute  the  operator  d/dt  =  -id),  we  obtain 


-icoy  =  ±ay 


(4.3.8) 


which  is  a  form  present  in  each  of  the  perturbation  equations.  Hence,  we  can  see  that  in 
the  electron  continuity  equation,  the  group 
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contains  both  a  growth  and  a  damping  term. 

Terms  which  contain  the  same  perturbed  property  as  the  unsteady  term  but 
additionally  are  multiplied  by  ±i  can  also  be  classified.  The  main  effect  of  these  terms  is 
on  the  frequency  of  oscillation,  and  are  hence  hereafter  referred  to  as  “frequency”  terms. 
This  can  be  clarified  by  considering  the  real  and  imaginary  parts  of,  for  example,  the 
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electron  continuity  perturbation  equation  separately,  and  examining  how  the  real  and 
imaginary  parts  of  co  are  affected  by  the  terms  in  the  equation  (multiplication  by  i  can 
also  be  thought  of  as  a  90  degree  shift  in  relative  phase).  For  example,  consider  the  first 
two  terms 


-icon,  +ikJJ ,  n„+  ... 


(4.3.10) 


of  the  continuity  perturbation  equation.  It  is  clear  in  this  example  that  the  second  term 
affects  only  the  real  component  of  co,  the  frequency  of  oscillation,  and  is  hence  a 
“frequency”  term. 

The  remaining  terms  which  cannot  be  classified  in  any  of  these  ways  have  in 
general  a  combined  effect  on  the  stability  and  dispersion  characteristics  of  the 
perturbations.  The  prime  factors  for  determining  the  effect  of  these  terms  on  the 
disturbances,  whether  stabilizing  or  destabilizing,  are  the  relative  phase  and  amplitude  of 
the  perturbations. 

Another  thing  to  note  is  that  in  the  momentum  equations,  the  coefficient  of  the 
unsteady  terms  (4/((l  +  for  ion  momentum  equation)  are  very  small,  signifying 

that  the  unsteady  terms  are  negligible  for  frequencies  much  smaller  than  the  ion-neutral 
collision  frequency.  Therefore,  the  unsteady  terms  in  the  momentum  equations  can  be 
ignored  for  most  frequencies  of  interest;  the  momentum  equations  would  then  describe  a 
balance  between  the  pressure,  field,  and  collisional  forces  due  to  the  perturbations.  This 
also  implies  that  for  these  frequencies,  the  stability  and  dispersion  characteristics  are 
mainly  determined  by  the  continuity  and  electron  energy  equations. 
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4.4.  The  Role  of  Ionization 

Examination  of  eqn.  4.2.6,  the  steady  state  electron  energy  equation,  reveals  that 
it  is  non-linear  in  Te  and  therefore  must  be  solved  for  a  given  electron  number  density 

and  current  density  by  using  an  iterative  technique.  Ignoring  for  now  the  recombination 
term  (the  last  term  within  the  brackets),  we  have  an  expression  for  the  rate  of  energy 
dissipation  for  the  given  current  and  electron  number  density.  Plotting  the  right  hand 
side  with  respect  to  Te ,  as  shown  in  Fig.  4.3,  reveals  a  maximum  value  of  energy  that  can 

be  dissipated  by  collisions  and  ionization,  which  approximately  corresponds  to  the 
maximum  rate  of  charged  particle  production.  This  maximum  occurs  at  the  condition: 

7:=Y  (4.4.1) 


Figure  4.3.  Variation  in  rate  of  energy  loss  due  to  ionization  and  coliisional  energy  exchange  with 
electron  temperature. 
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or,  dimensionally, 


(4.4.2) 


For  atomic  hydrogen  this  corresponds  to  a  value  of  52609  K,  or  4.53  ev. 

With  Fig.  4.3  we  may  begin  to  speculate  on  the  mechanisms  and  conditions  which 
might  cause  small  perturbations  to  become  unstable.  On  the  high  temperature  branch 
(temperatures  greater  than  that  corresponding  to  the  maximum  point),  a  small  increase  in 
(electron)  temperature,  by  way  of  a  perturbation,  results  in  a  decrease  in  the  rate  of 
energy  dissipation.  If  that  temperature  perturbation  is  accompanied  by  a  positive 
perturbation  to  the  energy  addition  rate  (from  the  joule  heating  and  pressure  work 
mechanisms),  then  the  result  may  be  an  initial  runaway  increase  in  the  temperature.  A 
runaway  situation  can  also  happen  if  the  opposite  process  occurs,  i.e  if  on  the  low 
temperature  branch,  a  decrease  in  the  electron  temperature  leads  to  an  increase  in  energy 
addition.  Another  possibility  is  that  an  increase  in  temperature  causes  an  increase  in 
energy  dissipation  which  is  smaller  than  the  increase  in  energy  addition.  The  slope  of  the 
curve  determines  this  increase  in  energy  dissipation,  which  implies  that  regions  of  the 
curve  with  shallower  slopes  would  tend  to  be  more  susceptible  to  instability.  Thus, 
although  the  point  of  inflection  may  seem  like  an  obvious  place  for  instability  to  occur, 
the  actual  location  of  the  demarcation  between  temporally  stable  and  unstable  situations 
is  unclear  and  must  be  determined  using  an  unsteady  analysis.  Note  that  this  simple 
thought  experiment  does  not  take  into  account  mechanisms  such  as  conduction  and 
pressure  work  which  also  have  an  effect  on  stability. 

The  ionization  rate  perturbation  is  given  in  the  right  hand  side  of  the  electron 
continuity  equation.  We  dimensionalize  this  term  by  multiplying  it  with  hR ,  the 
reference  production  rate.  The  resulting  dimensional  production  rate,  he ,  is  written  as 
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The  coefficient  of  this  term, 
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(4.4.3) 


(4.4.4) 


which  we  will  denote  the  “production  coefficient,”  provides  a  reasonably  simple  way  of 
quantifying  the  increase  in  charged  particle  production  due  to  perturbations  in  electron 
number  density  and  temperature.  For  higher  values  of  this  coefficient,  more  charged 
particles  are  produced  for  a  given  perturbation  in  electron  number  density  and 
temperature,  which  could  result  in  a  stronger  instability.  Care  must  be  taken  in  this 
interpretation,  however,  because  it  does  not  take  into  account  changes  in  relative  phase  of 
the  electron  number  density  and  temperature  perturbation  resulting  from  changes  in  k  . 
Also,  when  comparing  the  stability  of  different  steady  state  conditions,  it  must  be  noted 
that  the  terms  involving  the  perturbations  in  equations  4.3.1  -  4.3.5  also  depend  on  the 
steady  state  conditions. 

4.5.  Solution 

Eqns.  4.3.1  -  4.3.5  can  be  expressed  as  a  generalized  eigenvalue  problem,  with  a 
as  the  eigenvalue.  The  form  of  this  eigenvalue  problem  is 


|a  -  tniijy  =  0 


(4.5.1) 


where 
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y  = 


n. 


Ue  U;  E  Te  ,  the  perturbation  amplitudes 


and 


|a  -  a>i?J  is  the  coefficient  matrix  of  the  perturbation  eqns.  4.3. 1  -  4.3.5. 


The  solution  of  such  a  problem  can  be  handled  by  using  a  standard  routine  found  in 
many  commercial  matrix  analysis  programs  such  as  Matlab.  However,  it  is  also 
informative  to  consider  an  analytical  method: 

det[A-m\y=  0  (4.5.2) 

This  statement  can  then  be  reduced  to  an  algebraic  equation  using  a  symbolic  math 
program  such  as  Mathematica.  The  resulting  equation  is  a  third  order  polynomial,  which 
denotes  that  three  “eigenmodes”  are  possible  for  each  steady  state  condition.  Each 
eigenmode  has  a  distinct  value  for  d)  ,  the  complex  eigenvalue,  as  well  as  its  own  set  of 
perturbation  amplitudes. 

4.6.  Isothermal  waves 

The  solution  of  the  perturbation  equations  for  temporal  stability  reveals  that  there 
are  three  modes  which  can  propagate  and  grow  (or  decay)  in  this  initially  uniform 
plasma.  As  a  matter  of  convenience  we  wish  to  identify  these  modes  and  assign  a 
nomenclature  to  them  so  that  they  can  be  identified  for  different  steady  state  conditions. 
To  facilitate  this,  we  study  a  simplified  “isothermal”  model  which  gives  an  analytical 
solution  for  the  eigenvalues.  We  form  the  isothermal  model  from  the  set  of  perturbation 
equations  by  fixing  the  electron  temperature  at  the  temperature  of  the  heavy  particles, 
thus  removing  the  need  for  the  electron  energy  equation.  Setting  the  determinant  of  the 
resulting  system  to  zero,  we  obtain  the  following  quadratic  equation  for  (0: 
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which  we  identify  as  an  “ion”  disturbatnce,  since  the  phase  velocity  of  the  wave  ( coR/k ) 

is  very  close  to  the  ion  velocity.  We  shall  see  in  next  section  that  this  mode  is 
characterized  by  a  comparitively  large  ion  velocity  disturbance,  hence  additionally 
justifying  its  nomenclature. 

We  refer  to  the  other  mode  as  the  “pressure”  disturbance,  because  it  is  found  to  be 
more  closely  driven  by  gradients  in  the  electron  number  density  and  electron 
temperature.  Setting  the  inelastic  energy  loss  factor,  Sh,  to  a  very  high  value  and  solving 
the  full  set  of  perturbation  equations  also  reveals  these  two  isothermal  modes.  In  essence 
8h  affects  the  degree  of  departure  of  the  electrons  from  thermal  equilibrium.  When  the 
non-isothermal  case  is  investigated  (Sh  =  5,  a  typical  value  for  hydrogen)  as  described 
below,  a  third  mode  appears,  which  we  associate  with  the  electron  energy  equation  and 
thus  identify  as  a  “thermal”  wave.  This  assignment  is  also  related  to  the  fact  that  this 
third  mode  has  a  damping  rate  that  is  characterized  by  the  energy  dissipation  coefficients 


(4.6.3) 


in  the  electron  energy  perturbation  equation. 
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4.7.  Results 

4.7.1  Wave  Characteristics 

We  now  examine  the  dispersion  characteristics  of  the  three  wave  modes  identified 
from  the  analysis.  Of  these  three  modes,  the  thermal  wave  is  the  most  significant 
because  it  is  the  only  one  which  exhibits  unstable  behavior.  This  illustrates  the  point  that 
not  every  oscillatory  perturbation  behaves  in  the  same  way  with  regards  to  stability, 
reinforcing  the  notion  that  an  unsteady  analysis  must  be  used  to  determine  stability 
characteristics.  In  electrothermal  instability  theory,  a  positive  perturbation  in 
temperature  is  automatically  assumed  to  result  in  a  positive  perturbation  in  charged 
particle  number  density.  However,  we  shall  see  that  this  is  not  generally  true  and  that 
many  time  scale  and  spatial  factors  contribute  to  the  stability  characteristics  we  are 
interested  in. 

In  this  section  the  dispersion  characteristics  of  the  three  wave  modes  for  a  plasma 
at  6000  K  drawing  a  non-dimensional  current  density  of  /  =  400  are  presented.  This 
current  density  is  typical  of  the  near  anode  current  densities  calculated  in  low  power 
arcjet  thrusters.  The  remaining  parameters  are  again  given  in  Table  3.1.  A  set  of 
calculations  is  made  with  ne  =  0.5  and  Te  =3.13,  which  dimensionally  correspond  to  ~ 

8xl0‘8  m'3  and  18,000  K.  Most  information  that  is  needed  to  discuss  the  behavior  of  the 
three  modes  can  be  seen  in  the  plots  of  the  dispersion  characteristics,  growth/decay  rates, 
and  relative  phase  and  amplitude  vs.  wave  number,  which  are  presented  in  Figs.  4-13. 

We  begin  our  discussion  with  the  analysis  of  the  ion  wave  mode,  whose 
dispersion  and  growth/damping  rate  curves  are  given  in  Figs.  4.4  and  4.5,  respectively. 
By  examining  its  growth  rate,  we  see  that  it  remains  fairly  constant  at  a  value  of  about 
3xlO\  which,  when  dimensionalized  (divided  by  tR)  gives  8.4xl09  s'1,  the  ion-neutral 
collision  frequency.  Fig.  4.6,  which  shows  the  ion  velocity  amplitude  normalized  by  the 
electron  velocity  amplitude,  shows  that  the  ion  velocity  perturbation  amplitude  is  very 
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Figure  4.4.  Dispersion  curve  for  the  ion  mode.  Figure  4.5.  Growth  rate  curve  for  the  ion  mode. 


k 

Figure  4.6.  Ion/Electron  velocity  amplitude  for  the  ion  mode. 


nearly  equal  to  the  electron  velocity  perturbation  amplitude.  This  is  not  the  case  with  the 
other  two  wave  modes,  in  which  the  electron  velocity  perturbation  amplitude  is  much 
larger  than  that  of  the  ions.  These  two  observations  help  confirm  the  fact  that  the  “ion” 
wave  is  a  disturbance  that  is  in  fact  dominated  by  ion  motion.  Because  of  this,  the 
dominant  equation  for  this  mode  is  the  ion  momentum  equation,  which,  upon  removal  of 
field  and  pressure  force  terms  (since  they  are  small),  becomes  a  equation  relating  the 
unsteady  (inertia)  term  to  the  collisional  damping  term, 
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which  explains  the  damping  rate  and  behavior  of  this  mode. 

The  thermal  and  pressure  modes,  in  contrast,  are  modes  where  the  ion  velocity 
perturbation  amplitude  is  much  smaller  than  that  of  the  electron  velocity.  Thus,  the 
electron  number  density,  temperature,  and  velocity  perturbations  are  the  dominant 
properties  in  these  modes.  Of  the  two  modes,  the  thermal  wave  is  of  greater  interest 
because  of  its  potential  for  instability.  Figs.  4.7  -  4.10  show  the  dispersion  and 
growth/decay  characteristics  for  these  modes.  An  interesting  point  on  these  graphs  is  the 
k  =  0  point.  At  this  infinite  wavelength  condition,  many  of  the  terms  in  the  perturbation 
equations  drop  out,  leaving  only  certain  growth/damping  (such  as  the  energy  dissipation 
and  Joule  heating)  terms  to  balance  the  unsteady  terms.  In  addition,  these  growth/decay 
terms  remain  constant  for  all  values  of  k  ,  signifying  that  any  changes  in  character  for 
non-zero  values  of  k  are  due  to  the  remaining  terms,  which  include  the  mechanisms  of 
ion/electron  flux,  pressure  force,  energy  convection,  pressure  work,  and  electron  energy 
conduction. 

From  Figs.  4.7  -  4.10,  we  see  that  at  k  =  0,  the  frequency  of  oscillation  0)r  is  zero 
and  (x),  is  negative,  indicating  a  perturbation  which  decays  uniformly.  Therefore,  without 
the  influence  of  the  “gradient”  type  mechanisms,  the  uniform  underdense  plasma  is  still 
stable.  However,  as  seen  from  Fig.  10,  for  certain  values  of  k  ,  the  perturbation  is 
unstable.  The  instability  is  hence  dependent  on  these  mechanisms  for  instability, 
indicating  that  the  Joule  heating  perturbation  alone  is  not  sufficient  for  instability  in  these 
conditions. 

The  behavior  of  the  thermal  wave’s  stability  characteristics  can  be  partially  clarified  by 
examining  figure  1 1,  a  plot  of  the  relative  phase  between  the  electron  number  density  and 
temperature  perturbation.  For  the  thermal  wave,  we  see  that  the  relative  phase  starts  at 
180  degrees  for  k  =  0.  Eqn.  2.4.3  (the  continuity  pert  equation)  shows  that  the  term 
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Figure  4.7.  Dispersion  curve  for  the  pressure 
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Figure  4.8.  Growth  rate  curve  for  the  pressure 


k 

Figure  4.9.  Dispersion  curve  for  the  thermal 
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Figure  4.10.  Dispersion  curve  for  the 
thermal  mode.  The  dashed  line  denotes  a 
“critical”  condition  where  ne  =  0.9. 
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is  a  growth  term.  The  perturbation  to  charged  particle  production  due  to  the  electron 
temperature  perturbation  is 
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Figure  4.11.  Relative  phase  between 
the  electron  number  density  and 
temperature  for  the  thermal  mode. 


Figure  4.12.  Relative  amplitude 
between  the  electron  number  density 
and  temperature  for  the  thermal  mode. 
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Figure  4.13.  Relative  phase  between 
the  electron  number  density  and 
temperature  for  the  pressure  mode. 
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Figure  4.14.  Relative  amplitude 
between  the  electron  number  density 
and  temperature  for  the  pressure  mode. 


which  will  be  a  pure  damping  term  if  the  relative  phase  is  180  degrees,  and  a  pure  growth 
term  if  the  it  is  0  degrees.  At  90  or  270  degrees  the  temperature  perturbation  has  no 
effect  on  the  growth  or  decay  of  the  electron  number  density;  its  main  effect  is  on  the 
frequency  of  oscillation.  This  can  again  all  be  seen  by  considering  the  real  and 
imaginary  parts  of  the  electron  continuity  equation  separately,  and  examining  how  the 
real  and  imaginary  parts  of  CD  are  affected  by  the  different  terms  in  the  equation.  As 

k  increases,  the  relative  phase  decreases,  signifying  that  the  opposing  effect  of  the  two 
perturbations  on  the  production  perturbation  is  reduced  at  higher  wave  numbers. 
Therefore  the  growth  of  the  electron  number  density  is  due  to  the  increase  of  charged 
particles  caused  by  the  perturbation  in  the  electron  number  density.  We  thus  conclude 
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that  in  this  mode  the  temperature  perturbation  does  not  contribute  to  the  overall  increase 
in  charged  particle  production. 

The  relative  amplitude  of  the  number  density  perturbation  compared  to  the 
temperature  perturbation  also  increases  with  k  ,  as  shown  in  Fig.  4.12.  Again,  this  helps 
explain  why  the  thermal  mode  becomes  unstable  at  certain  k  . 

The  growth  rate  shows  a  maximum  value  at  k  ~  700.  This  maximum  exists 
because  of  the  influence  of  electron  energy  conduction,  which  is  clearly  a  damping 
mechanism  with  a  quadratic  dependence  on  the  wavenumber.  Hence,  at  high 
wavenumbers,  the  conduction  term  becomes  an  increasingly  important  stabilizing 
mechanism. 

The  same  phase  and  amplitude  plots  are  given  for  the  pressure  mode  in  Figs.  4.13 
and  4.14.  For  this  mode  we  see  that  the  relative  phase  between  the  number  density  and 
temperature  remains  for  the  most  part  fairly  close  to  180  degrees  with  changing  k  ,  with 
the  number  density  amplitude  also  dropping  with  respect  to  the  temperature  amplitude. 
This  behavior,  which  is  contrary  to  that  of  the  thermal  wave,  helps  explain  why  the  decay 
rate  grows  with  k  ,  and  why  the  mode  never  becomes  unstable. 

A  very  interesting  thing  to  point  out  is  that  if  the  pressure  work,  convection,  and 
conduction  terms  in  the  electron  energy  perturbation  equation  are  removed,  the  critical 
electron  temperature  increases  to  8.1  ev,  which  corresponds  to  a  dimensional  electron 
temperature  of  almost  50,000  K.  In  this  temperature  regime  Joule  heating  is  very  high 
and  the  slope  of  the  curve  in  Fig.  4.3  is  small.  This  signifies  that  Joule  heating  alone  is 
not  a  dominant  mechanism  in  the  instability  at  lower  electron  energies.  Because  the 
convection  term  acts  as  a  “frequency”  term,  the  main  contributer,  either  directly  or 
indirectly,  to  instability  must  be  the  pressure  work  perturbation,  where  its  main  influence 
may  be  its  effect  on  the  phase  relationships  between  the  perturbed  quantities.  It  should 
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be  remembered  that  the  pressure  work  term  does  not  destabilize  the  pressure  or  ion 
modes,  implying  that  it  can  also  act  as  a  damping  mechanism  in  these  modes. 

4.7.2.  Critical  conditions 

We  define  the  maximum  value  of  n.  (and  hence  minimum  T  )  where  the  curve 
of  a>;  for  the  thermal  mode  just  becomes  positive  at  any  value  of  k  as  the  “critical” 
value  of  ne*  .  Fig.  10  shows  the  growth/decay  rate  curves  for  two  values  of  he  :  the 

original  value  of  0.5,  and  the  critical  value  of  0.9.  Here  it  can  be  seen  that  as  the  number 
density  drops,  the  instability  becomes  stronger.  In  this  section  we  discuss  how  the 
critical  conditions  are  affected  by  varying  properties  of  the  background  gas  such  as  the 
ionization  energy,  pressure,  and  temperature. 

Examination  of  the  steady  state  equations  4.2.3  -  4.2.6  shows  that,  for  conditions 
where  recombination  can  be  neglected,  the  electron  velocity,  ion  velocity,  electric  field, 
and  electron  temperature  are  all  functions  of  the  parameter  J  /ne^ ,  s7 ,  9,  and  /i.  Under 

the  conditions  of  a  constant  background  temperature  Tx,  background  pressure  P,  gas  type, 
collision  frequencies  v,„  and  ven,  and  inelastic  energy  loss  factor  8h,  the  last  three 
parameters  are  constant.  For  this  idealized  condition  the  solution  to  the  steady  state 
equations  are  solely  determined  by  the  non-dimensional  ratio  of  current  density  to 
electron  number  density,  which  is  equal  to  the  non-dimensional  electric  field.  This  is  not 
quite  true  in  the  case  of  the  perturbation  equations  due  to  the  presence  of  ne  by  itself  in 

the  continuity  and  momentum  equations.  However,  we  have  found  that  despite  this  fact, 
the  critical  /  /n  value  is  essentially  constant  and  a  good  parameter  to  use  in 

determining  stability  characteristics.  Hence,  within  the  limits  of  this  simplified  model, 
smaller  current  densities  could  exhibit  instabilities  at  lower  electron  number  densities 
because  of  the  higher  rates  of  ionization,  while  larger  current  densities  can  become 
unstable  for  smaller  departures  from  Saha  equilibrium.  Fixing  this  non-dimensional 
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parameter  also  fixes  the  critical  electron  temperature,  which  can  similarly  be  viewed  as  a 
stability  criterion. 

4.7.2.I.  Ionization  Energy 

The  ionization  energy  would  seem  to  be  an  important  parameter  in  determining 
the  stability  characteristics  of  a  uniform  underdense  plasma  because  of  its  influence  on 
charged  particle  production  rates.  To  isolate  the  effect  of  the  ionization  energy  on 
stability,  stability  calculations  are  made  with  varying  values  of  £,  for  T  =  6000  K,  J  = 
1.1x10s  A/m2.  This  is  a  simple  way  learn  about  the  stability  characteristics  of  different 
gases,  though  variations  in  degenerate  states,  inelastic  collision  factors,  collision 
frequencies,  etc.  are  ignored  because  they  cloud  the  analysis  with  additional 
complications. 

The  effect  of  ionization  energy  on  stability  is  not  obvious  from  an  inspection  of 
the  perturbation  equations.  This  is  because  while  lower  ionization  energies  give  rise  to 
higher  production  rates  for  a  given  electron  number  density  and  temperature,  the  energy 
dissipated  by  each  charged  particle  pair  produced  is  lower.  Hence,  there  are  competing 
factors  whose  overall  influence  on  the  critical  conditions  is  not  clear.  The  critical 
conditions  are  calculated  and  presented  in  Figs.  4.15  and  4.16,  which  show  T.  and 

ne  versus  er  Here,  simple  relationships  are  seen  between  these  parameters.  Our  cutoff 

is  £j  =  8  ev,  below  which  the  critical  electron  temperature  drops  below  the  background 
temperature,  hence  entering  a  regime  where  there  is  steady  state  collisional  energy  gain. 
The  critical  electron  temperature  seems  to  vary  linearly  with  ep  with  higher  ionization 
energies  corresponding  to  higher  values  of  Te ,  .  The  critical  electron  number  density 

correspondingly  falls  with  increasing  ionization  energy.  The  size  of  variation  is  not 
large,  differing  by  roughly  a  factor  of  2  between  ionization  energies  ranging  from  8  to 
16.  Because  of  this,  the  ratio  J  /ne  ,  which  describes  the  critical  steady  state  electric 

field. 
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Figure  4.15.  Critical  non-dimensional 
electron  temperature  vs.  Ionization 

energy. 
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Figure  4. 1 6.  Critical  non-dimensional 
election  number  density  vs.  Ionization 
energy. 


increases  with  ionization,  implying  that  gases  with  higher  ionization  energies  tend  to  be 
more  stable. 


We  define  the  wavenumber  for  which  the  growth  rate  is  zero  as  the  critical 
wavenumber  kcnt .  In  performing  the  critical  conditions  analysis  we  find  that  the  critical 

wavenumber  stays  on  the  order  of  104  m'1  (with  a  tendency  to  decrease  with  increasing 
ionization  energy),  corresponding  to  a  critical  wavelength  of  0.6  millimeters. 

4.7.2.2.  Pressure  dependence  (T  =  6000,  £7=13.6  ev) 

The  changes  in  stability  characteristics  due  to  varying  background  (neutral 
species)  pressure  is  presented  in  this  section.  The  influence  of  changing  the  background 
pressure  is  felt  mainly  on  the  charged  particle  production  rate  and  the  charged  particle 
collision  frequencies,  both  of  which  depend  linearly  on  the  neutral  number  density.  The 
calculated  critical  conditions  are  presented  in  Table  4.1  for  background  pressures  ranging 
from  .1  atm  to  10  atm.  They  show,  interestingly,  that  the  pressure  has  no  effect  on  the 
critical  electron  temperature  and  number  density,  and  that  the  main  changes  are  in  the 
growth/damping  rates  and  critical  wavenumbers.  Increasing  the  background  pressure 
while  maintaining  a  fixed  background  temperature  is  equivalent  to  increasing  the  neutral 
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species  number  density,  which  results  in  increasing  the  number  of  particles  available  for 
ionization  and  collisions.  The  production  terms  therefore  scale  linearly  with  pressure,  as 
do  the  growth/damping  rates  (which  are  given  for  k=  0)  of  the  perturbations. 
Conductivity,  however,  scales  linearly  with  the  collision  frequencies,  so  that  for  a  fixed 
1  /ne^,  the  electric  field  also  scales  with  the  charged  particle  collision  frequencies. 

Hence,  at  lower  pressures,  critical  conditions  occur  at  lower  electric  fields  than  higher 
pressures.  However,  we  find  that  the  associated  growth/damping  rates  are  also  smaller, 


Pressure 

(atm) 

T  (kelvin) 

nt  rj  (m'3) 

kcnt  (m'1) 

Knt  (m) 

o 

11  p 

.1 

14500 

1.44  xlO19 

2770 

2.27  xl0“3 

5.61xl06 

.5 

14500 

1.44  xlO19 

13800 

4.55  xlO-4 

2.81xl07 

1 

14600 

1.43  xlO19 

27600 

2.28  xlO"4 

5.62xl07 

5 

14600 

1.42  xlO19 

138000 

4.55  xlO"5 

2.78  x10s 

10 

14500 

1.44  xlO19 

286000 

2.20  xlO-5 

Table  4.1.  Effects  of  varying  pressure  for  T  =  6000  K,  J  =  1  xlO5  A/m2,  e,  =  13.6  ev. 


signifying  that  the  instability  is  weaker.  For  higher  pressures,  then,  the  onset  of 
instability  may  be  delayed  a  small  amount,  but  the  instability  will  be  stronger  when  it 
occurs.  This  fact  may  be  more  useful  when  non-linear  effects  are  considered,  since  these 
effects  are  related  to  the  actual  magnitudes  of  the  disturbances. 


Another  parameter  which  varies  with  pressure  is  the  critical  wavenumber,  which 
increases  with  pressure.  This  means  that  the  critical  wavelength  decreases  with 
increasing  pressure,  which  may  also  be  an  important  factor  to  consider  when  taking  into 
account  physical  dimensions  of  plasma  devices,  as  well  as  the  extent  to  which 
electrostatic  effects  contribute.  The  explanation  behind  why  this  happens  has  to  do  with 
the  mean  free  path  decreasing  with  increasing  pressure,  which  affects  the  distance  over 
which  collisional  effects  are  important. 
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Along  with  the  changes  to  the  critical  wavenumber  we  find  that  the  frequencies 
associated  with  these  critical  conditions  also  vary  with  pressure.  Table  4.1  shows  that 
these  frequencies  increase  with  pressure,  beginning  from  Mhz  frequencies  at  P  =  .  1  atm 
to  100  Mhz  frequencies  at  10  atm.  This  observation  may  also  be  an  important 
consideration  in  finding  ways  to  keep  instabilities  from  forming. 

4.7.2.3.  Temperature  Dependence  (P  =  .5  atm,  e;=13.6  ev) 

We  can  speculate  about  the  effect  of  varying  heavy  species  temperature  by 
examining  the  general  conservation  equations,  eqns.  2.1.1  -  2.1.5.  Because  the  heavy 
species  energy  equations  are  neglected  in  this  study,  the  heavy  species  temperature  is  of 
secondary  interest.  The  heavy  species  temperature  nevertheless  does  have  an  influence  in 
each  equation  used  in  this  model.  In  the  continuity  equation,  this  temperature  affects  the 
neutral  species  number  density,  which  as  we  discussed  in  the  previous  section  affects  the 
generation  rate.  In  the  ion  momentum  equation,  its  presence  is  felt  in  the  pressure 
gradient  term.  In  the  electron  energy  equation,  the  difference  between  electron  and  bulk 
temperatures  helps  determine  the  collisional  energy  exchange  between  the  two  types  of 
particles.  And  finally,  the  electron-neutral  and  ion-neutral  momentum  exchange 
frequencies,  which  affect  both  the  momentum  and  electron  energy  equations,  depend  on 
the  neutral  species  number  density,  which  again  is  related  to  both  the  temperature  and 
pressure. 

The  influence  of  the  bulk  temperature  can  be  seen  in  Table  4.2,  which  describes 
its  effect  in  terms  of  the  critical  electron  temperature  and  electron  number  density  for  a 
fixed  current  density  J  =  1.1x10s.  The  results  show  that  the  background  temperature  does 
have  an  effect  on  the  critical  conditions;  the  critical  electron  temperature  and  number 
density  both  increase  slightly  with  bulk  temperature,  though  the  changes  are  not  great. 
Because  of  this,  \J  /n.  )  ,  and  hence  Ecrit  also  increase  with  bulk  temperature,  which 

can  be  explained  by  the  fact  that  a  stronger  field  is  required  to  heat  the  electrons  to  the 
critical  electron  temperature. 
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T  (kelvin) 

Te  (kelvin) 

nc  (m"3) 

Ecrit  (N/C) 

Krit  (m) 

4000 

14300 

1.35  xlO19 

1.09  xlO8 

3.05  xlO"1 

5000 

14400 

1.39  xlO19 

9.51  xlO7 

6000 

14500 

1.45  xlO19 

8.37  xlO7 

4.33  xlO"4 

7000 

14600 

7.40x10’ 

4.76  xlO-4 

Table4.2.  Effects  of  varying  background  temperature  for  P  =  .5  atm,  J  =  lxlO5  A/m2,  e,  =  13.6  ev. 


The  bulk  temperature  also  has  an  effect  on  the  wavenumber  at  which  the  instability  first 
appears.  With  increasing  temperature,  this  critical  wavenumber  decreases,  which 
indicates  that  the  critical  wavelength  is  increasing.  The  frequency  of  these  perturbations 
also  change  slightly,  while  remaining  at  about  the  same  order  of  magnitude  (in  the  Mhz 
range). 


The  analysis  in  this  chapter  has  provided  us  with  a  simplified  way  to  gain  insight 
into  the  possible  mechanisms  responsible  for  instability.  By  examining  this  uniform 
situation  first  we  have  examined  frequencies,  wavelengths,  energy  equation  dependence, 
the  importance  of  ionization,  and  introduced  the  role  of  phase  and  amplitude 
relationships  which  are  also  necessary  to  understand  a  more  realistic  stability  analysis,  as 
we  shall  see  in  the  following  chapter.  In  addition,  insight  into  how  background 
conditions  affect  the  results  of  the  analysis  has  been  gained,  which  again  will  enhance 
our  approach  to  interpreting  the  results  of  the  non-uniform  problem. 
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Chapter  5. 

Stability  Analysis  -  Non-Uniform  Plasma  Case 


We  now  remove  the  artificial  constraint  of  a  uniform  plasma  and  analyze  a  more 
physically  realistic  situation  which  takes  into  account  the  spatially  varying  background 
plasma  in  the  electric  boundary  layer.  To  accomplish  this,  we  recall  the  analysis  which 
was  discussed  in  chapter  3,  where  the  steady  state  profiles  for  the  electric  boundary  layer 
were  calculated.  As  outlined  in  Chapter  2,  these  conditions  are  analyzed  for  stability  by 
modeling  the  time-dependent  response  of  the  boundary  layer  to  small  planar 
perturbations  which  have  the  form  of  Eqn.  2.4.2,  which  upon  substitution  results  in  the 
perturbation  equations  2.3.2  -  2.3.6. 

5.1.  Additional  simplifications 

In  the  actual  solution  of  the  non-uniform  stability  problem,  it  is  very  desirable  to 
keep  the  number  of  equations  as  small  as  possible.  One  reduction  can  be  made  by 
recalling  the  far-field  boundary  condition  that  all  perturbations  die  out  far  from  the 
anode,  and  combining  this  boundary  condition  with  eqn.  2.3.3  ( V  -  J  -  0,  ignoring  two- 
dimensional  effects).  This  leads  to  the  simple  result 

t=t  (5.1.1) 


which  states  that  the  perturbation  in  the  overall  current  density  perturbation,  J ,  is  not 
only  constant  but  equal  to  zero  throughout  the  entire  boundary  layer.  Rearranging  eqn. 
5.1  results  in  the  statement  that 
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which  can  then  be  substituted  into  the  original  set  to  eliminate  Ur  The  other 

simplification  that  can  be  made  is  by  adding  the  electron  and  ion  momentum  perturbation 
equations  (eqns.  2.3.4  and  2.3.5)  to  remove  the  electric  field  perturbation.  The  final 
remaining  set  thus  consists  of  three  equations  with  three  unknowns;  the  original  electron 
continuity  perturbation  equation  (eqn.  2.3.2),  the  electron  energy  perturbation  equation 
(eqn.  2.3.6),  and  the  following  statement  for  the  electron  velocity  perturbation  (ignoring 
electron  inertia): 
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5.2.  Solution  method 


(5.1.3) 


The  actual  solution  of  the  eigenvalue  problem  can  be  obtained  numerically  by 
three  different  methods:  finite  differences,  the  shooting  method,  or  a  spectral  type 
method.  Simple  finite  differences  are  used  in  this  study  to  discretize  the  equations  and 
generate  the  matrix  from  which  the  eigenvalues  are  calculated.  The  advantages  of  using 
the  finite  difference  method  are  mainly  in  simplicity  and  ease  of  programming.  The 
disadvantages  are  slow  convergence  and  the  return  of  a  large  number  of  “spurious”  roots, 
which  are  characteristic  of  the  discretization  matrix  only  and  have  no  relation  to  the 
physical  solution.  The  number  of  these  spurious  roots  grows  with  the  size  of  the  matrix, 
which  is  a  good  motivator  for  reducing  the  set  of  equations  to  the  smallest  number 
possible.  Distinguishing  the  physical  roots  from  the  spurious  roots  requires  the 
calculation  of  eigenvalues  from  matrices  of  different  sizes  (which  represent  different  grid 
spacings),  and  finding  the  eigenvalues  which  are  invariant  with  grid  size  (these  are  the 
“physical”  roots). 

Discretization  of  the  domain  can  be  done  either  with  a  stretched  or  uniform  grid. 
In  this  study,  both  were  used  to  help  confirm  the  existence  of  the  eigenvalues.  Results 
are  given  here  only  for  uniform  grid  calculations. 
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The  discretized  form  of  the  three  perturbation  equations  are  written  for  each  grid 
point  and  form  a  generalized  eigenvalue  problem  with  a  coefficient  matrix  and 

eigenfunction  vector,  as  discussed  in  chapter  2.  The  A  (coefficient)  and  B  matrices, 
however,  are  now  square  matrices  with  3N  rows  and  columns,  where  N  is  the  number  of 
grid  points.  We  note  that  this  will  result  in  much  larger  matrices  than  those  encountered 
in  the  uniform  analysis.  Numerical  solution  times  of  generalized  eigenvalue  problems 
scale  approximately  with  N3,  which  again  reinforces  the  need  for  keeping  the  number  of 
equations  and  grid  points  small.  In  general,  the  maximum  number  of  eigenmodes 
calculated  is  3N,  the  size  of  the  matrix;  in  practice,  the  actual  number  of  eigenmodes  is 
slightly  below  this  value. 

Theoretically,  the  physical  roots  should  converge  to  their  values  as  the  number  of 
grid  points,  and  hence  accuracy,  is  increased.  In  this  study,  the  cases  are  calculated  with 
a  certain  number  of  grid  points,  and  then  re-calculated  with  twice  as  many  grid  points  as 
before.  Physical  eigenvalues  should  be  revealed  for  both  cases  and  be  fairly  equal, 
showing  a  pattern  of  convergance  with  decreasing  mesh  size.  Spurious  eigenvalues  will 
not  display  this  type  of  behavior,  and  also  will  tend  to  have  the  additional  characteristic 
that  their  eigenfunctions  display  strong  spatial  variations  which  bear  a  resemblance  to  so- 
called  “noise”  plots. 

The  three  perturbation  equations  represent  a  fourth  order  set  of  linear  ordinary 
equations,  which  require  four  boundary  conditions  for  solution.  As  discussed  in  Chapter 

A  A 

2,  we  require  that  ne ,  Ue ,  and  Te  dissapear  in  the  far-field.  The  remaining  boundary 
equation  is  the  enforcement  of  eqn  2.6.1,  the  recombination  boundary  condition  at  the 
anode.  The  perturbation  form  of  this  boundary  condition,  obtained  by  the  usual  method 
of  substituting  the  sum  of  steady  state  and  perturbation  quantities,  is 
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at  the  anode  wall.  We  note  that  it  is  not  necessary  to  include  the  questionable  condition 
of  a  zero  gradient  for  the  electron  temperature  perturbation  at  the  anode  for  this  non- 
uniform  analysis. 

5.3.  Results 

5.3.1  High  current  case,  unstable  mode 

The  conditions  which  we  examine  for  stability  are  again  those  given  in  Table  3.1, 
supporting  a  non-dimensional  current  density  of  J  =  40,  which  corresponds  to  a 
dimensional  current  density  of  J  =  1.1  xlO4  A/m2.  The  boundary  layer  profiles  for  these 
conditions  are  also  presented  in  Chapter  3.  The  current  density  chosen  was  identified  in 
Eskin’s  thesis  as  an  approximate  transition  point  from  the  diffuse  mode  to  the  constricted 
mode,  though  it  is  not  stated  where  this  value  originates  from.  At  higher  current 
densities,  the  eigenfunction  gradients  are  sharper  and  hence  a  greater  spatial  resolution  is 
needed  to  accurately  resolve  them. 

Performing  the  calculations  using  a  uniform  grid  with  100  grid  points,  over  a 
domain  from  J  =  0  to  x  =  .03,  we  find  that  the  resulting  eigenvalues  span  over  a  wide 
range  of  values.  A  significant  number  of  these  eigenvalues  are  not  realistic  because  they 
have  (dimensional)  real  components  which  are  smaller  than  10'10  seconds,  which  is  on  the 
order  of  a  100  year  period!  These  roots  are  clearly  spurious  and  can  immediately  be 
discarded.  The  remaining  set  of  eigenvalues,  which  still  contain  some  spurious  modes, 
are  displayed  on  a  log-log  scale  of  their  imaginary  vs.  real  parts  in  Fig.  5.1. 

We  first  note  that  the  Fig.  5.1  is  symmetric  about  the  C0r  axis.  The  eigenvalues 
which  differ  only  by  a  sign  in  the  real  part  actually  represent  the  same  mode.  Another 
thing  to  note  is  that  there  are  quite  a  few  roots  where  ar  lies  in  the  106-108  range  for 

which  fit  is  positive.  A  discussion  of  these  roots  will  follow;  we  will  first  focus  our 
attention  on  the  lone  root  where  3r  lies  in  the  104  range  and  has  a  positive  imaginary 
component.  This  root  has  a  calculated  eigenvalue  of  fi)  =  3.75xl04  +L32xl04/  for  100 
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grid  points.  Repeating  the  calculation  with  200  grid  points  reveals  that  the  root  takes  on 
a  value  of  3.85  x  104  + 1.26  x  104  i ,  and  with  300  grid  points  (which  takes  about  a  full  day 
to  run  on  a  Sun  Sparc  20  workstation),  m  =  3.91  xlO4 +1.23x10 Ai,  which  indicates 
fairly 


log  10  (Or 

Figure  5.1.  Log-log  plot  of  the  real  vs.  the  imaginary  parts  of  the  eigenvalues  for  100  grid  points. 

good  convergance.  We  will  take  this  as  an  indication  that  the  root  is  not  spurious  and 
hence  physical.  It  should  be  noted  that  the  dimensional  value  of  the  frequency 
(i cor  =  a>r/ tj j)  for  this  unstable  eigenvalue  is  around  106  s'1,  or  about  1  Mhz. 

Examination  of  the  eigenfunctions  of  this  unstable  mode  provides  insight  on  how 
the  perturbations  vary  spatially.  The  amplitude  of  the  eigenfunctions  of  the  three 
parameters  which  are  directly  calculated  -  ne,  Uc,  and  Tc  -  are  given  in  Fig.  5.2  for  the  200 
grid  point  calculation,  with  every  other  data  point  (100  points  total)  plotted  to  reduce 
clutter.  Each  eigenfunction  amplitude  is  normalized  with  its  maximum  value  so  that  the 
three  may  be  plotted  on  the  same  graph  for  comparison.  We  find  that  the  electron 
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velocity  and  temperature  amplitudes  show  a  smooth,  almost  exponential,  decay  towards 
zero,  with  the  velocity  exhibiting  a  shorter  scale  length  than  the  temperature.  The 
number  density  amplitude  shows  a  maximum  near  the  anode  before  its  smooth  decay, 
which  may  be  a  result  of  the  finite  rate  recombination  boundary  condition  (eqn.  5.2.1). 
The  eigenfunction  amplitudes  of  the  properties  which  are  “backed  out”  after  the  main 

problem  is  solved  -  U, f ,  E  -  as  well  as  illustrative  properties  such  as  the  electron/ion 

flux,  F j ,  and  flux  gradient  dre  /  dx  perturbations  -  are  shown  in  Figs.  5.3  and  5.4, 

respectively.  The  relative  phase  plots  of  all  of  these  eigenfunctions  are  given  together  in 
Fig.  5.5. 

We  see  from  these  plots  that  near  the  anode,  the  eigenfunctions  for  the  ion 
velocity  and  electron  flux  are  not  very  smooth,  and  are  characterized  by  steep  gradients 
which  appear  as  discontinuities  because  of  the  grid  scale  used.  This  is  again  due  to  the 
finite  rate  recombination  boundary  condition  imposed  at  the  anode.  Written  in  another 
way,  it  imposes  that 


(-2k,n„-U,jK=n„V,  (5.3.5) 

at  x  =0.  This,  along  with  eqn.  5.1.1,  are  statements  that  attempt  to  model  the  effect  of 
the  anode  while  maintaning  the  one-dimensional  and  quasineutral  simplifications. 
However,  this  is  physically  unrealistic  because  of  the  existence  of  a  sheath  in  front  of  the 
anode.  When  electrostatic  and  two-dimensional  effects  are  taken  into  account,  there  is  no 
longer  the  requirement  to  set  the  electron  and  ion  flux  perturbations  equal  to  each  other. 
Also,  the  use  of  a  single  recombination  coefficient  to  describe  the  highly  complex  anode- 
plasma  interaction  is  a  gross  simplification.  In  order  to  more  realistically  treat  the 
problem,  the  ion  and  electron  flux  should  be  decoupled,  and  separate  reaction  rate 
coefficients,  each  with  their  own  dependencies  on  temperature  and  number  density, 
should  be  used  to  describe  specific  reactions,  in  addition  to  taking  two-dimensional  and 
electrostatic  effects  into  account.  However,  while  acknowledging  the  shortcomings  of 


Figure  5.2.  Normalized  Amplitudes  for  the  eigenfunctions  ne ( □),  Ue(»),  Te  (+).  The  maximum  values 
of  the  amplitudes  are  1.88xl07,  0.353,  and  3.52x10  *,  respectively. 


Figure  5.3.  Normalized  Amplitudes  for  the  eigenfunctions  ne( □),  J/;(v),  E(  0).  The  maximum 
values  of  the  amplitudes  are  1.88xl0'7, 1.33x10  ',  and  8.91x10",  respectively. 
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eqn.  5.2.4,  we  will  still  rely  on  it  to  model  the  anode  boundary  condition  in  a 
mathematically  simple  way. 

The  discontinous  behavior  of  the  ion  velocity  and  flux  near  the  anode  disappears 
for  low  values  of  ksj.  This  can  be  seen  in  eqn.  5.3.5,  where  higher  values  of  ks  force  the 
ion  velocity  and  electron  number  density  to  be  180  degrees  out  of  phase.  In  the  limit  of 
setting  ka  to  zero,  the  ion  flux  perturbation  is  forced  into  phase  with  the  number  density 
perturbation  (because  the  steady  ion  velocity  at  the  anode  is  negative  for  non-saturated 
cases),  which  gives  the  ion  velocity  eigenfunction  a  smoother  variation  and  removes  the 
steep  gradient  behavior  in  the  ion  velocity  amplitude.  In  the  limit  of  the  one-dimensional 
quasineutral  model,  it  seems  that  the  ion  velocity  and  flux  are  better  behaved  near  the 
anode  for  low  values  of  ksi,  which  physically  corresponds  to  lower  steady  state  ion  fluxes 
at  the  anode  (closer  to  saturation  conditions). 

At  a  small  distance  away  from  the  anode,  we  find  that  every  one  of  the 
eigenfunctions  are  well  behaved.  The  electron  velocity  disturbance  amplitude  is  many 
orders  of  magnitude  greater  than  that  of  the  ions,  which  as  we  saw  in  Chapter  4  indicates 
that  this  is  not  an  “ion’  type  disturbance.  The  phase  of  the  electron  number  density  and 
velocity  are  very  nearly  (but  not  exactly)  equal.  If  they  were  exactly  equal,  the  ion 
velocity  perturbation  would  be  either  in  phase  or  180  degrees  out  of  phase  with  these  two 
properties  (as  seen  by  eqn.  5.1.2),  which  is  not  the  case.  The  phase  of  the  electron 
number  density  perturbation  also  does  not  change  very  rapidly,  making  a  90  degree 
change  from  the  anode  to  approximately  x  =  .02,  where  its  amplitude  has  almost  died 
out.  This  shows  that  the  number  density  perturbation  does  not  exhibit  a  strong  spatial 
oscillation,  and  that  the  perturbation  is  mostly  monotonic.  This  observation  may  be 
important  when  considering  the  non-linear  processes  of  arc  transition.  We  will  see  that 
for  another  set  of  conditions  (the  low  current  density  case)  the  number  density  phase  has 
a  stronger  spatial  variance,  which  essentially  divides  the  disturbance  into  regions  where 
the  electron  number  density  perturbations  are  negative  and  positive. 
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A  physical  explanation  for  the  instability  exhibited  by  this  mode  can  be  deduced 
by  examining  the  nature  of  the  steady  state  boundary  layer.  By  considering  electron 
particle  conservation,  we  recall  that  the  steady  state  boundary  layer  is  in  a  state  of 
charged  particle  production,  which  is  balanced  out  by  the  flux  gradients  which  carry  the 
produced  particles  towards  the  anode,  where  they  recombine.  By  first  assuming  for  the 

moment  that  the  steady  state  production  rate  is  fixed  (that  is,  h  =  0),  it  is  possible  to 
hypothesize  that  a  perturbation  to  the  overall  flux  of  charged  particles,  perhaps  due  to  a 
fluctuation  in  the  electric  field,  can  upset  the  particle  balance  and  result  in  a  build-up  of 
charged  particles  in  the  region.  This  process  could  then  conceivably  lead  to  a  local 
increase  in  conductivity,  which  has  been  speculated  as  the  transition  initiating  mechanism 
in  electrothermal  instability  theory. 

This  flux  instability  mechanism  can  also  be  aided  by  h .  For  the  case  under 
consideration,  a  perturbation  to  the  overall  electron/ion  flux  results  in  a  positive 
perturbation  to  the  charged  particle  number  density,  which  grows  due  to  a  combination 
of  the  flux  disturbance  and  increase  in  particle  production.  Important  physical  insight 
can  be  gained  by  examining  eqn.  2.4.3  (the  continuity  perturbation  equation).  From  this 
equation  we  conclude  that  the  perturbation  flux  gradient  is  a  major  factor  in  determining 
the  growth  or  decay  characteristics  of  the  disturbance.  By  comparing  the  eigenfunctions 
of  the  electron  number  density  and  flux  gradient  (which  is  calculated  by  a  simple 
difference  formula),  we  can  determine  the  role  that  the  flux  gradient  is  playing  in  the 
instability.  Fig.  5.4  shows  that  the  two  eigenfunctions  start  out  relatively  in  phase,  and 
actually  come  into  phase  at  x  ~  .0025.  In  this  region,  the  flux  perturbation  gradient  is 
acting  to  damp  out  the  disturbance  by  moving  particles  out  of  this  near  anode  region. 
However,  there  is  a  contribution  to  the  growth  of  the  number  density  perturbation  from 

h  ,  which  maintains  the  instabililty  in  this  region.  The  “production  coefficient,”  which 
was  introduced  in  the  previous  chapter,  is  fairly  high  (about  104)  here  and  thus  indicates 
that  many  additional  particles  are  being  rapidly  produced  in  the  region  due  to 
perturbations  in  electron  density  and  temperature.  Further  away  from  the  anode,  where 
the  production  coefficient  is  smaller,  the  two  eigenfunctions  go  out  of  phase  by  about 
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120  degrees,  indicating  that  the  flux  gradient  begins  to  aid  the  instabilty.  A  relative 
phase  of  180  degrees  would  indicate  that  the  flux  gradient  is  purely  growth  term,  though 
anything  from  90  to  270  degrees  indicates  a  net  influx  of  particles  and  hence  a  partial 
role  as  a  growth  mechanism. 

This  characterization  of  the  instability  as  a  “flux’"  instability  is  something  that  was 
not  exhibited  in  the  previous  “quasi-uniform”  analysis.  We  recall  in  that  analysis  that  the 
steady  state  flux  gradients  were  ignored,  which  prevented  the  flux  perturbation  from 
becoming  a  major  factor  in  determining  stability  characteristics.  Instability  was  related 
to  the  electron  energy  equation  as  a  disruption  to  the  energy  balance  of  Joule  heating, 
pressure  work  and  conduction,  which  are  more  closely  tied  with  electrothermal  instability 
theory.  In  the  present  non-uniform  analysis,  the  electron  energy  equation  does  not 
appear  to  be  the  deciding  factor  in  determining  the  instability  of  the  boundary  layer. 
Hence,  it  seems  that  net  ionization  rate,  rather  than  current  density,  is  the  important 
factor  in  current  stability.  This  is  reinforced  by  the  following  case,  which  is  a  stability 
analysis  of  a  very  low  current  boundary  layer  in  the  state  of  Joule  cooling. 

5.3.2.  Low  current  case,  unstable  mode 

The  steady  state  profiles  for  a  boundary  layer  with  the  same  conditions  as  the 
previous  case,  with  the  exception  that  J  -  A  (about  110  A/m2),  are  also  given  in  figs. 
3.1  -  3.8.  For  this  case  we  find  that  the  steady  state  electron  temperature  falls  as  the 
anode  is  approached,  indicating  a  negative  electric  field  and  the  boundary  layer  being  in 
a  state  of  “Joule  cooling”.  Despite  this  fact,  the  boundary  layer  is  still  in  a  state  of 
ionization,  although  the  production  rates  are  much  lower  than  in  the  previous  case. 
Again  we  can  imagine  a  flux  perturbation  in  this  boundary  layer  which  can  cause  the 
number  density  of  the  electrons  and  ions  to  grow,  though  we  would  expect  the  growth  to 
be  weaker.  The  stability  analysis  reveals  the  existence  of  an  unstable  mode  with 
ft)  =  1.9  x  104  +  4.4  x  10 3  i ,  which  has  a  growth  rate  about  a  factor  of  3  smaller  than  the  J 
=  40  case.  Figs.  5.6  -  5.9  show  the  plots  of  the  eigenfunction  amplitudes  and  phase 
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relationships,  in  a  format  similar  to  those  given  for  the  higher  current  case.  We  find  that 
while  there  are  similarities,  there  are  also  some  distinguishing  characteristics  for  this  low 
current  mode  which  could  potentially  be  important  when  considering  the  dynamics  of 
anode  spot  formation. 

It  is  obvious  by  considering  the  steady  state  temperature  profile  that  the 
production  coefficient  is  always  going  to  be  very  close  to  unity.  Hence,  the  production 

perturbation  h  is  not  going  to  have  a  strong  influence  on  any  instability.  If  a  mode  is 
going  to  be  unstable,  it  should  be  entirely  due  to  the  flux  perturbation.  This  fact  is 
confirmed  by  examining  the  phase  plot  for  this  case,  which  shows  that  the  relative  phase 
between  the  electron  number  density  and  flux  gradient  perturbations  is  always  greater 
than  90  degrees. 

Another  interesting  difference  for  this  case  is  that  the  phase  of  the  number  density 
perturbation  changes  very  rapidly  near  the  anode,  making  a  180  degree  change  from  the 
anode  at  x  =  .01.  This  is  in  contrast  to  the  higher  current  case,  where  the  number 
density  phase  did  not  vary  significantly  past  90  degrees.  This  indicates  that,  for  this  low 
current  case,  the  number  density  perturbation  has  more  of  a  spatially  wave-like  character, 
where  regions  are  divided  into  positive  and  negative  perturbations.  Physically,  this 
signifies  that  charged  particles  are  initially  only  being  shuffled  around  by  the 
disturbance,  and  that  there  are  very  few  new  particles  being  produced.  These 
characteristics  may  have  important  implications  when  considering  transition  dynamics, 
though  any  further  conclusions  would  require  verification  from  more  complex  models. 

We  note  that  for  this  low  current  case,  we  do  not  see  the  discontinous  behavior  of 
the  ion  velocity  and  electron  flux  which  was  found  in  the  previous  case.  This  is  due  to 
the  fact  that  for  a  large  portion  of  the  boundary  layer,  the  steady  state  ion  velocity  is 
negative,  resulting  in  a  more  continuous  variation  for  the  eigenfunctions  in  question. 


Stability  Analysis  -  Non-uniform  Plasma  Case 


ca 

0 

- 

A 

o  A 

□  o 

□  A 
°D  \ 

Figure  5.6.  Normalized  Amplitudes  for  the  eigenfunctions  ne  ( □ ),  Ue  (  0  ),  Te  (A).  The  maximum 
values  of  the  amplitudes  are  1.09xl0'\  0.471,  and  1.72xl0'3,  respectively. 


Figure  5.7.  Normalized  Amplitudes  for  the  eigenfunctions  ne(U\  dre/dx  (•),  JTe{ A).  The  maximum 
values  of  the  amplitudes  are  1.09xl0'5,  0.21,  and  6.11xl0\  respectively. 
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The  results  of  the  low  current  case  seem  to  indicate  that,  for  this  particular  mode, 
the  electric  boundary  layer  instability  is  based  at  least  partially  on  the  upsetting  of  the 
balance  of  steady  state  charged  particle  production  and  flux.  Joule  heating  seems  to  have 
a  very  minor,  if  any,  role  in  the  stability  characteristics.  This  fact  is  confirmed  by  the 
fact  that  raising  the  current  density  by  an  order  of  magnitude  also  has  a  relatively  minor 
effect  on  the  growth  rate.  Specifically,  a  stability  analysis  for  the  J  =  400  case  shows 
that  co  =  4.0  x  1 04  + 1.38  x  104 i  for  this  mode,  which  does  not  indicate  a  significant 
increase  in  growth  rate  from  the  7  =  40  case.  Part  of  the  reason  for  this  is  that  the  peak 
steady  state  electron  temperature  (which  occurs  at  the  anode)  is  1.83,  as  opposed  to  1.76 
for  the  7  =  40  case.  Hence,  the  steady  state  production  rate  h  and  the  production 
coefficient  do  not  increase  very  significantly  from  increasing  the  current  density  by  an 
order  of  magnitude.  For  boundary  layers  which  do  have  a  strong  current  density  to 
ionization  rate  relationship,  however,  we  might  expect  a  stronger  correlation  between  the 
current  density  and  the  growth  rate  to  exist. 

5.3.3.  Other  modes 

As  mentioned  previously,  the  stability  analysis  of  all  current  densities  discussed 
so  far  turned  up  a  number  of  eigenvalues  which  also  had  positive  imaginary  values.  It  is 
difficult  to  say  whether  these  modes,  which  have  non-dimensional  frequencies  and 
growth  rates  on  the  order  of  107,  are  physical.  This  is  because  the  eigenfunctions  show 
very  sharp  gradients  throughout  the  entire  region,  with  many  modes  exhibiting  a  very 
“noisy”  pattern.  The  number  of  grid  points  needed  to  properly  resolve  such  gradients  is 
large  and  are  very  numerically  expensive  to  calculate.  As  a  result,  to  keep  within  the 
scope  of  this  study,  we  will  not  focus  any  more  attention  on  these  modes.  Their  steep 
gradients  and  non-convergance  give  serious  doubts  to  their  existence,  and  more  elaborate 
methods  are  needed  to  analyze  them  to  determine  if  they  are  indeed  physical. 

We  now  turn  our  attention  to  the  examination  of  the  stable  modes.  Moving  up  in 
frequency  from  the  physical  unstable  mode  discussed  earlier,  there  exist  a  set  of 
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eigenfunctions  with  non-dimensional  frequencies  and  damping  rates  in  the  105  range, 
about  an  order  of  magnitude  greater  than  the  unstable  mode.  As  an  illustrative  example, 
we  choose  to  examine  a  mode  where  co  ~  1.65  xlO5  -3.29  xlO5/,  which  has  its 
eigenfunction  amplitudes  and  phases  plotted  in  Figs.  5.10  -  5.13.  From  these  plots  it  can 
be  seen  that  at  these  higher  frequencies,  some  of  the  eigenfunctions  display  a  more 
spatially  oscillating  character.  This  is  illustrated  by  the  changes  and  periodic  variations 
of  the  perturbation  phase.  A  perfectly  straight  phase  line  would  correspond  to  an  exact 
sinusoidal  variation  in  space.  The  amplitude  plots  seem  to  indicate  a  “beat”  behavior  in 
some  of  the  eigenfunctions.  An  interesting  thing  to  note  is  that  for  a  certain  region  of  the 
boundary  layer  (3c  >  .0075  ),  the  flux  gradient  perturbation  is  almost  exactly  in  phase 
with  the  number  density  perturbation.  This  helps  explain  why  the  mode  is  stable,  though 
we  suspect  that  the  spatially  oscillating  character  of  the  eigenfunction  also  plays  a  strong 
role.  Again,  it  should  be  noted  that  none  of  the  eigenfunctions  suffer  from  any 
discontinuities. 

The  remaining  type  of  eigenmode  that  can  be  resolved  occurs  for  eigenvalues 
with  damping  rates  very  close  to  3x  108 ,  which  is  a  strong  indication  that  these  are  “ion” 
type  disturbances  first  encountered  in  the  uniform  analysis.  Examination  of  the 
eigenfunctions  for  these  modes  again  strengthen  the  notion  that  the  higher  frequency 
modes  tend  to  be  strongly  oscillating  “noise”  plots  which  we  cannot  comment  any  further 
on.  Lower  frequency  ion  modes,  however,  are  very  well  behaved.  An  example  of  this  is 
the  mode  corresponding  to  ty  *  7.05 xlO5  -3.0 xlO8/ ,  which  is  plotted  in  Figs.  5.14  - 
5.17.  A  comparison  of  the  ion  and  electron  velocity  magnitudes  shows  that  their 
amplitudes  are  very  nearly  equal,  confirming  the  notion  that  this  disturbance  is 
dominated  by  ion  momentum.  Fig.  5.17,  which  displays  the  relative  phase  plot,  clearly 
indicates  by  the  linear  variation  that  all  of  the  eigenfunctions  are  sinusoidal.  This  figure 
also  shows  that  the  number  density  and  flux  gradient  perturbations  are  almost  exactly  in 
phase,  which  helps  explain  the  stability  properties  of  the  mode. 


Figure5.13.  Relative  phase  (in  degrees)  for  the  eigenfunctions  ne(U),  Ue  (o),  Te  (+),  I/,(v),  E  (0), 
dTe  /  dx  (  a  ).  Note  that  the  electron  density  and  velocity  perturbations  are  approximately  in  phase,  and 
that  phase  of  Te  is  not  plotted  because  it  has  the  same  phase  as  the  ion  velocity  perturbation. 


Figure  5.14.  Normalized  Amplitudes  for  the  eigenfunctions  ne(\ □),  Ue  (°  ),  Te  (A),  ion  mode.  The 
maximum  values  of  the  amplitudes  are  1.09xl0's,  0.471,  and  1.72xl0'3,  respectively. 


Figure  5.15.  Normalized  Amplitudes  for  the  eigenfunctions  ne(p)7  dFe/dx{ o),  re( A),  ion  mode.  The 
maximum  values  of  the  amplitudes  are  1.09x10  s,  0.21,  and  6.11xl0'4,  respectively. 
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Figure  5.16.  Normalized  Amplitudes  for  the  eigenfunctions  Ui  (□),  E  (°),  ion  mode.  The  maximum 
values  of  the  amplitudes  are  1.55x1  O'2  and  .131,  respectively. 


Figure  5.17.  Relative  phase  (in  degrees)  for  the  eigenfunctions  ne  (□ ),  Ue  (  °),  Te (+),  U{  ( v ),  E  ( 0  ), 
dTe  /  dx  (  a  ),  Te  (*),  ion  mode.  Note  that  the  electron  density  and  velocity  perturbations  are 
approximately  in  phase. 
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It  is  interesting  to  note  that  among  the  different  modes,  the  electron  temperature 
perturbation  maintains  a  fairly  consistent  quality.  This  could  possibly  be  interpreted  as 
additional  evidence  that  the  electron  energy  equation  does  not  have  a  large  role  to  play  in 
the  stability  characteristics  of  the  electric  boundary  layer. 

5.4.  Model  Limitations 

Because  of  the  short  length  scales  which  characterize  the  perturbations,  we  must 
be  aware  of  the  quasineutral  limitations  of  our  model.  We  have  already  discussed  the 
possible  effects  of  the  sheath  on  the  anode  boundary  conditions  in  a  previous  section. 
We  must  also  consider  whether  Poisson’s  equation  (eqn.  2.1.6),  and  hence 
quasineutrality,  is  being  violated  by  the  perturbations.  The  condition  for  satisfying 
Poisson’s  equation  for  the  perturbations  is: 


o  ClJZi 

ze  —  «  zne  (5.4.1) 

ax 

where  the  parameter  z  is  used  to  illustrate  that  the  perturbation  vector  y  represents  a 
linear  “basis”  vector  for  the  domain  of  allowable  perturbations.  However,  it  is  clear  that 
z  drops  out  and  has  no  effect  in  this  analysis,  which  allows  the  usage  of  any  convenient 
perturbation  value  to  determine  the  applicability  of  the  quasineutral  assumption. 

As  a  rough  approximation,  we  can  use  the  perturbation  magnitudes  to  give  an 
order  of  magnitude  “test”  for  Poisson’s  equation.  For  our  current  case,  we  find  that  e 

has  a  value  of  4.6  x  10”5 ,  and  that  the  scale  length  of  the  electric  field  perturbation  can  be 
very  conservatively  estimated  to  be  0.001.  Substituting  the  maximum  values  for  the 
electric  field  and  electron  number  density  perturbations,  we  find  that  the  left  hand  side  of 
eqn  5.4.1  is  about  1%  the  value  of  the  right  hand  side,  which  indicates  that 
quasineutrality  seems  to  be  reasonably  satisfied  for  these  perturbations. 
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We  must  always  be  aware  that  our  assumption  of  a  one-dimensional  problem  also 
introduces  approximations.  However,  a  more  rigorous  two-dimensional  treatment  brings 
a  whole  host  of  theoretical  and  numerical  complications.  To  help  understand  the 
approximations  being  made  with  the  1-D  assumption,  we  examine  the  differential  form 
of  Ampere’s  Law: 


V  x  B  =  |i0 


J  +£0 


dE' 

dt  j 


(5.4.2) 


Taking  the  divergance  of  both  sides  gives 


V-7  =  -e0-(v.£) 


(5.4.3) 


Expanding  this  into  2-D  form  by  assuming  radial  symmetry,  we  get 


1  d  ,  x  dJx  d 

~r~dSY  r^  +  lk~lk 


1  d  f  _  v  dEx 


(5.4.4) 


In  laboratory  situations,  we  expect  that  in  general,  real  perturbations  will  have  multi¬ 
dimensional  components.  However,  for  our  idealized  model  of  perfectly  planar 
perturbations,  we  can  neglect  the  radial  terms.  Hence,  after  non-dimensionlization  and 
applying  perturbations,  we  are  left  with 


(1+tf  at,  g 

4fl  tRoR  dx 


~  2x10  9  ico 


(5.4.5) 


for  the  conditions  stated  in  Table  3.1.  Using  an  order  of  magnitude  estimate  and 
assuming  that  the  gradients  on  both  sides  of  the  equation  have  the  same  scale  length  of 
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variation,  we  can  get  an  estimate  of  the  frequency  for  which  there  exists  a  net  current 
perturbation  (when  the  right  hand  side  of  the  equation  is  not  negligible)  as 


co  -5x10' 

max 


max 

E 

max 

~107 


(5.4.6) 


which  is  satisfied  for  the  unstable  mode  described  in  section  5.3.1. 

5.5.  The  role  of  the  boundary  layer  character  on  stability 

It  is  apparent  from  the  analysis  presented  that  the  model  lacks  the  ability  to  define 
a  boundary,  or  demarcation  point,  between  conditions  of  stable  and  unstable  behavior  for 
the  situation  described.  Differences  in  mode  structures  can  be  examined  to  speculate 
their  effects  on  current  mode  transition;  however,  the  task  of  testing  such  speculation 
requires  the  inclusion  of  more  mechanisms  and  realistic  conditions.  Although  the 
process  of  unraveling  the  mystery  of  arc  constriction  will  take  much  additional  study, 
additional  insight  and  observations  can  still  be  made  based  on  previous  studies  and  their 
steady  solutions  to  help  construct  future  research  efforts. 

As  discussed  in  the  introduction,  most  theoretical  treatments  have  concentrated  on 
steady  state  analysis  to  provide  insight  into  the  subject  of  arc  mode  transition.  These 
studies  have  shown  that  steady  current  transfer  in  plasmas  tends  to  take  place  along 
certain  pathways  which  minimize  energy  input  (the  Steenbeck  Minimum  Principle, 
which  stipulates  essentially  that  a  current  takes  the  path  of  least  resistance)  and  that  the 
boundary  layer  voltage  drop  reaches  a  maximum  at  a  certain  current  density  (the  point 
where  the  differential  impedance  goes  to  zero). 

For  the  results  presented  in  this  chapter,  we  again  compare  boundary  layers 
supporting  non-dimensional  current  densities  of  0.4  and  40,  respectively.  In  the  lower 
current  case  we  pointed  out  that  the  region  was  in  a  state  of  Joule  cooling.  In  Eskin’s 
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thesis,  such  regions  were  shown  to  lie  along  the  “positive  impedance”  section  of  the 
current-voltage  characteristic.  This  basically  means  that  in  this  current  regime, 
increasing  the  current  density  would  require  an  increase  in  the  voltage  supplied  to  the 
boundary  layer.  However,  in  the  negative  impedance  current  regime,  where  the 
boundary  layer  is  in  the  state  of  Joule  heating,  the  opposite  is  true  due  to  the  large 
numbers  of  charged  particles  being  produced.  This  can  simply  be  thought  of  as  the  result 
of  an  overall  increase  in  conductivity;  because  there  are  more  charged  particles  to  carry 
the  current,  the  voltage  needed  to  transfer  the  current  decreases.  This  line  of  reasoning 
cannot  be  used  in  describing  the  regime  of  positive  impedance,  however,  because  of  the 
fact  that  the  concept  of  “conductivity”  really  only  applies  to  a  simple  uniform  plasma 
where  the  current  and  field  both  act  in  the  same  direction.  For  the  Joule  cooling  cases, 
the  current  and  field  flow  in  opposite  directions  from  each  other. 

It  is  possible  that  these  distinctions  between  the  low  and  high  current  cases  will 
cause  planar  perturbations  to  interact  differently  with  the  steady  boundary  layer.  From 
Figs.  5.5  and  5.9,  we  find  that  the  electric  field  perturbation  is  about  180  degrees  out  of 
phase  with  the  number  density.  For  the  higher  current  case,  J  =40,  where  the  boundary 
layer  is  characterized  by  Joule  heating,  this  means  that  the  increase  in  the  electron 
number  density  due  to  the  perturbation  is  coupled  with  a  decrease  in  the  electric  field. 
This  corresponds  to  an  increase  in  conductivity  and  decrease  in  the  potential  in  the  region 
which,  from  the  above  reasoning,  could  be  interpreted  as  a  possible  transition  mechanism 
from  a  diffuse  attachment,  higher  voltage  mode  to  a  constricted  attachment,  lower 
voltage  mode.  For  the  lower  current  case,  the  increase  in  electron  number  density  again 
corresponds  to  a  decrease  in  electric  field.  However,  for  this  case,  the  background 
electric  field  is  negative,  which  signifies  that  the  overall  field  would  become  more 
negative.  The  potential  of  the  layer  would  hence  increase  as  the  charged  particles  are 
pushed  with  more  force  away  from  the  anode.  Heuristically,  it  seems  unlikely  that  such 
perturbations  would  result  in  arc  constriction  under  these  conditions. 
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The  above  speculation  can  only  be  confirmed  or  disproved  with  the  use  of  a  more 
general  model  which  takes  into  account  nonlinearity,  electrostatic  fields,  and  two- 
dimensionality.  In  the  present  analysis  there  is  no  way  to  analyze  the  behavior  of  the 
perturbations  beyond  their  initial  growth  phase. 

5.6  The  role  of  convection  on  stability 

As  discussed  in  the  introduction,  bulk  fluid  convection  has  been  experimentally 
shown  to  act  as  a  stabilizing  mechanism  [12].  As  a  result,  convection  may  play  a 
significant  role  in  the  stabilization  of  the  arc  in  arcjet  thrusters,  as  well  as  other  plasma 
devices.  In  the  previous  steady  state  analyses,  it  was  shown  that  convection  reduces  the 
ionization  rate  requirements  in  the  electric  boundary  layer  by  contributing  an  additional 
flux  of  charged  particles  [10]  and  hence  has  an  impact  on  the  electric  field  and  voltage 
characteristics.  These  observations  may  help  to  explain  the  stabilizing  effect  of  bulk 
fluid  motion,  since  it  is  hypothesized  in  this  study  that  the  strength  of  the  instability 
depends  on  the  net  ionization  rate  and  electrical  character  of  the  boundary  layer. 
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Chapter  6.  Conclusions  and  Recommendations 


The  model  presented  in  this  thesis  represents  a  fundamental  approach  to 
identifying  and  understanding  some  of  the  principal  mechanisms  responsible  for  arc 
mode  transition  and  breakdown.  The  results  of  the  model  provide  initial  insight  into  this 
process  as  well  as  a  solid  foundation  to  base  more  detailed  analyses  on.  This  chapter 
summarizes  the  results  of  the  study,  relates  the  findings  to  the  previous  studies  done  on 
the  subject,  and  provides  suggestions  for  future  studies  on  arc  stability. 

6.1.  Summary 

Our  stated  objective  for  this  study  was  to  model  the  time  dependent  response  of 
small  oscillating  disturbances  in  the  near-anode  region  and  to  identify  mechanisms  which 
affect  the  growth  or  decay  of  these  disturbances.  To  this  end  we  began  with  the 
calculation  of  the  steady  state  characteristics  of  the  near- anode  region  for  a  current 
carrying  plasma,  an  approach  found  in  many  more  recent  studies.  The  set  of  conditions 
we  chose  to  model  were  based  on  the  approximate  conditions  (based  on  previous 
numerical  simulations  of  arcjets)  found  in  the  near-anode  current  attachment  region  in 
arcjet  thrusters.  Many  simplifications,  such  as  treating  transport  and  reaction  coefficients 
as  constant  and  ignoring  secondary  species,  were  made.  Hypothesizing  that  mode 
transition  is  a  transient  phenomenon,  we  formulated  a  stability  analysis  based  on  small 
perturbations  in  the  plasma  properties  with  a  sinusoidal  time  dependence.  Assumptions 
were  again  made  to  simplify  the  stability  analysis,  the  most  important  one  being  the 
modeling  of  the  anode  and  perturbations  with  a  planar  geometry. 

To  better  obtain  an  understanding  of  the  results  we  anticipated  with  our  analysis, 
we  paused  briefly  to  perform  an  analysis  based  on  an  unphysical  assumption  -  that  of  a 
uniform  underdense  plasma.  By  using  this  simplification  we  again  significantly  reduced 
the  mathematical  complexity  of  the  analysis,  while  gaining  important  insight  into  the 
process  of  interpretating  the  results.  More  specifically,  we  were  able  to  more  simply 
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familiarize  ourselves  with  the  concepts  of  phase  and  amplitude  and  how  they  helped 
explain  the  stability  characteristics  of  the  calculated  modes.  We  were  also  able  to  begin 
to  identify  the  roles  of  many  mechanisms  such  as  Joule  heating,  ionization,  and 
conduction  on  the  growth  or  decay  rates  of  the  perturbations.  Classifications  were  given 
to  the  modes  which  resulted  from  the  analysis;  an  example  of  this  was  the  “ion”  mode, 
which  we  found  to  be  a  stable  mode  where  ion  motion  dominated  the  perturbation 
characteristics.  Finally,  we  were  able  to  gain  additional  insight  into  how  basic  properties 
such  as  bulk  temperature,  pressure,  and  the  gas  type  could  influence  stability. 

The  final  analysis  of  this  thesis  was  that  of  the  non-uniform  case,  which  took  into 
account  the  effect  of  the  gradients  inherent  to  the  electric  boundary  layer  on  arc  stability. 
From  this  analysis  we  found  that  the  model  was  not  able  to  define  a  specific  demarcation 
point  between  stable  and  unstable  conditions.  However,  the  influence  of  specific 
mechanisms  on  the  stability  characteristics  of  the  model  became  more  clear.  It  was  seen 
for  a  boundary  layer  characterized  by  Joule  heating,  significant  ionization  and  elevated 
electron  temperature  that  an  unstable  mode  existed  which  was  best  described  by  runaway 
charged  particle  production.  An  unstable  mode  was  also  found  in  a  very  low  current 
boundary  layer  case,  though  upon  closer  examination  important  differences  in  character 
from  the  higher  current  case  were  found.  This  low  current  instability  was  described  as  a 
“flux”  instability,  where  growth  in  the  charged  particle  density  was  the  result  of  a  flux 
imbalance  which  redistributed  the  existing  charged  particles,  rather  than  a  significant 
production  of  new  particles.  The  effects  of  these  differences  on  the  actual  process  of  arc 
transition  was  speculated  on  and  the  verification  left  to  future  study.  In  addition,  modes 
which  were  stable  were  also  identified  and  examined  to  help  explain  their  stability 
characteristics. 

6.2.  Comparisons  to  previous  studies 

We  recall  that  many  studies  which  involved  vacuum  arc  stability  found  that  anode 
vaporization  was  a  significant  factor  in  mode  transition.  It  was  hypothesized  that  this 
vaporization  resulted  in  a  significant  number  of  ionizable  neutral  particles,  which  then 
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provided  a  starting  point  for  the  transition  to  an  anode  spot.  Because  of  this,  factors 
which  are  important  to  anode  vaporization  -  anode  material,  geometry  and  heat  flux  into 
the  anode  -  were  identified  as  having  a  strong  influence  on  critical  current  densities.  In 
the  higher  pressure  arcs  we  are  interested  in,  the  number  density  of  the  neutral  species  is 
assumed  to  be  fixed  parameter.  However,  we  found  in  our  model  that  significant 
ionization  is  also  an  important  mechanism  for  arc  instability.  It  is  also  interesting  to  note 
that  the  voltage,  number  density,  and  temperature  fluctuations  we  would  expect  to  see  in 
the  higher  pressure  cases  have  been  observed  in  vacuum  arc  transition,  implying  that 
similarities  between  the  two  pressure  regimes  exist. 

In  our  review  of  experimental  high  pressure  arc  studies,  we  found  a  common 
thread  that  bulk  convection  had  a  significant  role  in  maintaining  a  stable,  diffuse  arc. 
Meeks  predicted  in  her  numerical  model  a  transition  between  diffusion-dominated  and 
convection-dominated  behavior  at  a  certain  bulk  fluid  flow  velocities.  She  found  that 
operating  a  diffuse  arc  under  convection-dominated  conditions  reduced  the  tendency  of 
the  plasma  to  ionize  in  the  boundary  layer  region,  which  we  have  found  to  change  the 
character  of  the  perturbations.  Because  the  model  presented  here  neglects  convection,  we 
believe  that  it  captures  and  reflects  much  of  the  important  qualitative  behavior  of  real 
diffuse  arcs  in  diffusion-dominated  conditions. 

6.3.  Future  directions  for  study 

We  believe  that  the  phenomenon  of  anode  spot  formation  is  inherently  multi¬ 
dimensional;  therefore  the  addition  of  at  least  a  radial  dimension  to  the  model  would 
enhance  its  ability  to  more  accurately  capture  the  transition  process.  In  addition,  non¬ 
linear  effects  will  surely  become  important  as  the  amplitude  of  the  perturbations  grow, 
and  should  investigated  in  future  studies.  To  better  understand  the  anode-plasma 
interaction’s  effect  on  stability,  more  realistic  boundary  conditions  which  take  into 
account  the  various  reaction  mechanisms  and  phenomena  occuring  at  the  anode  would 
provide  improved  coupling  between  the  plasma  and  the  electrode.  Specifically, 
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including  the  sheath  would  clarify  the  behavior  of  the  perturbations  in  non-neutral 
regions. 

The  inclusion  of  convection  into  the  steady  state  profiles  would  better  relate  the 
modeling  work  to  the  previously  discussed  experimental  observations.  A  parametric 
study  on  the  effects  of  varying  background  pressure,  temperature,  and  gas  types  would 
provide  better  understanding  of  the  current  model’s  capabilities.  And  finally,  the 
possible  effects  of  an  externally  applied  or  internally  generated  magnetic  field,  which 
some  researchers  believe  to  have  a  stabilizing  influence,  could  be  incorporated  into  a 
newer  multi-dimensional  model. 
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